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We introduce an event-by-event perturbative-QCD -|- saturation + hydro ("EKRT") framework 
for ultrarelativistic heavy-ion collisions, where we compute the produced fluctuating QCD-matter 
energy densities from next-to-leading order perturbative QCD using a saturation conjecture to 
control soft particle production, and describe the space-time evolution of the QCD matter with 
dissipative fluid dynamics, event by event. We perform a simultaneous comparison of the central¬ 
ity dependence of hadronic multiplicities, transverse momentum spectra, and flow coefficients of 
the azimuth-angle asymmetries, against the LHC and RHIC measurements. We compare also the 
computed event-by-event probability distributions of relative fluctuations of elliptic flow, and event- 
plane angle correlations, with the experimental data from Pb+Pb collisions at the LHC. We show 
how such a systematic multi-energy and multi-observable analysis tests the initial state calculation 
and the applicability region of hydrodynamics, and in particular how it constrains the tempera¬ 
ture dependence of the shear viscosity-to-entropy ratio of QCD matter in its different phases in a 
remarkably consistent manner. 


I. INTRODUCTION 

The main goal of ultrarelativistic heavy-ion collisions 
at the Large Hadron Collider (LHC) and the Relativistic 
Heavy-Ion Collider (RHIC) is to understand collectivity 
in the strong interaction sector of the Standard Model, 
and determine the properties such as temperature de¬ 
pendences of the shear and bulk viscosities in the differ¬ 
ent phases of QCD matter. Currently, with an increas¬ 
ing number of heavy-ion bulk observables from the LHC 
and RHIC to investigate, and with significant theoreti¬ 
cal developments over the last decade both in computing 
the produced initial state from QCD and in describing 
the subsequent space-time evolution with dissipative fluid 
dynamics event by event, one is now more concretely ap¬ 
proaching this ambitious goal. 

Bulk (1 ow-pt) observables - hadronic multiplicities, 
transverse momentum (px) spectra and especially the 
Fourier coefficients (u„) of their azimuth-angle distri¬ 
butions - measured in heavy-ion collisions at the LHC 
and RHIC, offer compelling evidence of a formation 
of a strongly collective locally nearly-thermalized low- 
viscosity hot QCD matter which undergoes both the 
quark-gluon plasma (QGP) and hadron resonance gas 
phases. For recent reviews, see mi- The measurements 
are remarkably consistent with describing the space-time 
evolution of the formed system with dissipative relativis¬ 
tic fluid dynamics Consequently, relativistic fluid 

dynamics has established its role as a cornerstone in the 
analysis of heavy-ion bulk observables. 

One of the clearest signals of a collective behavior of 
the matter produced in nuclear collisions is the emergence 
of azimuthal asymmetries of the hadron transverse mo¬ 


mentum spectra. In the fluid-dynamical limit the spatial 
inhomogeneities of the initial state are translated by the 
pressure gradients into the momentum space anisotropies 
of the spectra, and the effectiveness of this transition is 
essentially determined by the properties of the matter 
itself. It has turned out that the shear viscosity of the 
QCD matter strongly affects the final observed asymme¬ 
tries, and therefore the measured azimuthal structure of 
the transverse momentum spectra (quantified by the 
coefficients) gives the most direct constraints to the shear 
viscosity. 

As external input for solving the fluid-dynamical equa¬ 
tions of motion, one needs to know the QCD equa¬ 
tion of state (EoS) as well as event-by-event fluctuat¬ 
ing initial conditions for the spatial distribution of en¬ 
ergy (or entropy) density, the initial flow of the matter, 
and the starting time (space-time surface) for the evo¬ 
lution. The observable final-state momentum distribu¬ 
tions of hadrons are obtained by computing the hadronic 
momentum distributions at the decoupling of the sys¬ 
tem and accounting for resonance decays after that. To 
model the dynamics of hadron gas, including its dissipa¬ 
tion, decoupling and also the resonance decays, the fluid- 
dynamical evolution may also be coupled to a hadron 
cascade simulation at a suitably chosen space-time hy¬ 
persurface. Such hybrid approaches have been developed 
e.g. in [HIIISHII], see Ref. for a review. Common 
to the different dissipative fluid-dynamical settings, how¬ 
ever, is that the initial conditions play a crucial role in 
determining the uncertainties to the QCD matter prop¬ 
erties like the shear viscosity. 

A traditionally used way to get a hold on the initial 
conditions (see e.g. m m m 11313) is to assume the 
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initial energy (or entropy) densities to be a function of the 
Glauber model binary-collision and/or wounded-nucleon 
transverse densities and exploit the measured centrality 
dependence of various bulk observables (and more de¬ 
tailed observables such as relative EbyE fluctuations of 
Vn) for fixing the initial conditions in different centrality 
classes. A drawback in this is that there is essentially no 
predictability in the initial conditions when moving from 
one collision energy to another but the data fitting must 
be done for each cms-energy separately. Without consid¬ 
ering the QCD dynamics responsible for the initial gluon 
and quark production one does not have enough dynam¬ 
ical control over the formation time of the hot system, 
either. In this case, the freedom in re-iterating the initial 
conditions complicates the determination of the matter 
properties such as the temperature dependence of the 
shear viscosity. 

The uncertainties in the initial conditions, and thereby 
also in the QCD-matter viscosity determination, can be 
reduced if instead of fitting one can compute the initial 
conditions in a QCD-based framework. Steps into this 
direction include, e.g., the following approaches: 

• In the "IP-Glasma" initial conditions [3 m], one 
combines the impact parameter dependent color-glass- 
condensate (CGC) saturation model (=IP-Sat model) 
with a pre-thermal classical evolution of the glasma gluon 
fields. Combined with the MUSIC fluid-dynamics code 
El [13, such initial conditions have been particularly suc¬ 
cessful in explaining, e.g., the relative EbyE fluctuations 
of Vn measured by ATLAS [H] and ALICE EH). This ap¬ 
proach reproduces the measured Vn and u„(pr) system- 
atics very well with an effective constant shear-viscosity- 
to-entropy ratio 77 /s = 0.12 at RHIC and 0.2 at the LHC 

HI- 

• The Monte Carlo version of the Kharzeev-Levin- 
Nardi ("MC-KLN") model [SDHll], which is based on 
the CGC and kr factorization but where no pre-thermal 
evolution of the produced gluons is considered, has been 
used for obtaining the initial conditions in, e.g., |H1 [53] 
for the VISHNU hybrid code [TOl [31] . This setup gives 
a very good description of the measured multiplicities, 
Pt spectra and elliptic flow of bulk hadrons at RHIC 
and LHC assuming a constant viscosity-to-entropy ratio 
in the QGP, 77 /s = 0.16 [33]. As discussed in [S], com¬ 
paring the RHIC results obtained with the MC Glauber 
and MC-KLN initial conditions, one has arrived at an 
uncertainty interval 1 < 47r(77/s)QGP < 2.5. 

• The perturbative QCD + saturation model, often 
referred to as the Eskola-Kajantie-Ruuskanen-Tuominen 
("EKRT") model [3S], whose EbyE Next-to-leading or¬ 
der (NLO) extension we introduce here, combines the 
idea of the dominance of multiple few-GeV partonic jets, 
minijets, in high energy nuclear collisions |3S1[37] with 
a conjecture of saturation of gluon production to sup¬ 
press the non-perturbative particle production |58| . The 
original EKRT model |35l [39] , where the NLO effects in 
minijet transverse energy production snisi] were only 
partially accounted for, and where only ideal 1 D and 


1+1 D Bjorken hydrodynamics was applied, predicted 
the charged hadron multiplicities surprisingly correctly 
for central collisions both at the LHC [32] and RHIC 
|43| . Also the pt spectra of identified bulk hadrons at 
RHIC were reproduced very well [441 US] - For predic¬ 
tions of elliptic flow in this framework, using 2+1 D ideal 
fluid dynamics, see [23] for RHIC and [35] for the LHC. 

It is worth recalling here that the centrality depen¬ 
dence of multiplicities predicted by the EKRT model m 
was first thought not to agree with the RHIC measure¬ 
ments, see e.g. |481 135] . However, an excellent match 
with the data was eventually realized when the same (op¬ 
tical) Glauber model was used to calculate the number of 
participants also in the data analysis EOJEI] — compare 
Fig. 23(a) in and Fig. 22 (left) in with Fig. 4 in 
m- This observation also motivated us to develop the 
model further. In [52| we verified, albeit still using ideal 
hydrodynamics and leading order (LO) minijet cross sec¬ 
tions, that the EKRT model was able to reproduce well 
the bulk (1ow-pt) part of the LHC charged hadron pT 
spectrum in central Pb+Pb collisions. In [S3] the model 
was then consistently brought to NLO, its model param¬ 
eters were more precisely specified, the parameter cor¬ 
relations and propagation of nuclear parton distribution 
function (nPDF) uncertainties [S3] into the final multi¬ 
plicities were studied, and the predictive power of the 
model was demonstrated. 

Viscous fluid dynamics in the context of the NLO- 
improved EKRT model was introduced in |55| . where 
we performed a simultaneous analysis of the centrality 
dependence of charged hadron multiplicities, px spec¬ 
tra and elliptic flow, simultaneously for Pb+Pb colli¬ 
sions at the LHC and Au+Au at RHIC. The consis¬ 
tency of the EKRT results with the experimental data 
suggested, in terms of a linear parametrization assuming 
a minimum of 77 /s at T = 180 MeV, that 0.12 < 77 /s < 
0.12 + (0.18/320)(r/MeV- 180) in the QGP phase, and 
77 /s(T) = 0.12 - (0.20/80)(T/MeV - 180) in the hadron 
gas phase. Even though such a general behavior, a rising 
slope in T in the QGP is expected on the basis of lattice 
QCD [35] and a decreasing one in the hadron gas on the 
basis of kinetic theory 123. we also had to conclude in m 
that an equally good overall fit to the studied RHIC and 
LHC data can be obtained with a constant 77 /s « 0.20. In 
magnitude, this agrees with earlier studies [3 HT21 [551162| . 

To pin down the possible temperature dependence of 
77 /s in the different phases of QCD matter, further con¬ 
straints from analysing more detailed observables are 
needed. With this goal in mind, and especially for ac¬ 
cessing higher Fourier flow-coefficients and their EbyE 
analysis, we introduce here for the first time an EbyE 
framework to the NLO-improved pQCD + saturation + 
viscous fluid dynamics model [55] • The following issues 
and observables are considered in what follows: 

In Sec. H we define the 2+1 D equations of motion of 
longitudinally boost-invariant dissipative Israel-Stewart 
type transient fluid dynamics we use in this study, specify 
the parameters in our fluid dynamical setup, and discuss 
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the applicability of fluid dynamics in general. We also 
specify the Sf corrections to the local eqnilibrium parti¬ 
cle momentnm distribntion fnnctions, which are applied 
in the compntation of final state particle momentum dis¬ 
tributions at deconpling. Unfortunately, we are not yet 
capable of performing a full statistical global analysis of 
the LHC and RHIC heavy-ion measnrements to extract 
rj/s{T) and its uncertainty limits. However, as a step 
towards such an analysis, in order to demonstrate how 
sensitive (or, in some cases insensitive) the considered 
LHC and RHIC observables are to the shear viscosity, we 
study here the set of different parametrizations of rj/s{T) 
given in Sec. |HC[ 

In Sec. HI we explain in detail how the NLO-improved 
pQCD + satnration initial conditions are obtained EbyE, 
first addressing the infrared (IR) and collinear (CL) safe 
NLO calcnlation of minijet transverse energy and the 
conjectnre of saturation to obtain the satnration momen¬ 
tum psat locally in each transverse location. Accounting 
for the geometrical fluctuations of nucleon positions and 
exploiting the exclusive electroproduction measurement 
of J/'0 mesons at HERA [S3], we bnild np the initial 
gluon clonds in the colliding nuclei. The key point en¬ 
abling the EbyE framework in onr case in practice, is the 
scaling of psat with the prodnct of nnclear thickness func¬ 
tions of the colliding nuclei |551 IM] . From the local psat 
we then form the EbyE EKRT initial conditions, i.e., 
the energy densities and formation times locally in the 
transverse plane, addressing also the "pre-thermal" evo¬ 
lution to a constant longitudinal proper time tq = 0.2 fm 
at which we start the fluid dynamical simulation. Cen¬ 
trality selection and entropy prodnction during the fluid- 
dynamical evolntion in the EbyE case are demonstrated. 
Examples of the EKRT initial energy densities and ec¬ 
centricities vs. centrality are given, and the effects of the 
key parameters in our framework on the centrality de¬ 
pendence of the initial state entropy, eccentricities, and 
Psat are charted. 

Section IV snmmarizes the definitions of the flow- 
related observables, the Vn coefficients from 2-, 3- and 
4-particle cumulants, and event-plane angle correlations, 
which we compnte in the EbyE EKRT framework and 
compare with experimental data. 

Section V contains the resnlts from the new EbyE 
EKRT framework. We perform a systematic mnltiob- 
servable analysis, simultaneonsly for Pb+Pb collisions 
at the LHC and for the An+Au collisions at the top- 
energy of RHIC. We stndy the centrality dependence of 
charged hadron mnltiplicities, pt spectra, average pt’s of 
the identified bulk hadrons, and in particnlar the charged 
hadron flow coefficients and event-plane angle correla¬ 
tions. Also the probability distributions of the relative 
flnctuations of elliptic flow (<5^2) are computed and com¬ 
pared with LHC data as well as with the relative initial 
eccentricity fluctuations {5e2,5ei^2) in our EbyE EKRT 
setup. The necessity of fluid dynamics in understanding 
the centrality systematics of these quantities is demon¬ 
strated. 


In Sec. VI we discuss the applicability limits of the 
pQCD + saturation + fluid dynamics framework in the 
light of the computed flow coefficients and event-plane 
angle correlations, demonstrating the effects of the Sf 
corrections and showing where these effects start to be¬ 
come too large to be trnsted. 

The main conclusions from our new EbyE EKRT 
framework, discnssed in Sec. VH, can be snmmarized as 
follows: The computed centrality dependence of charged 
hadron mnltiplicities, low-py spectra, flow coefficients at 
the LHC and RHIC, and even the event-plane angle cor¬ 
relations at the LHC all agree very well with experimen¬ 
tal data for rj/s{T) = paraml, i.e. when rj/s(T) is mod¬ 
estly rising with T in the QGP and where r]/s(T) re¬ 
mains small in the hadron gas phase, see Fig. An 
equally good overall agreement is obtained with a con¬ 
stant r]js = 0.2. In particular, we strongly emphasize the 
necessity for a simultaneous analysis of LHC and RHIC 
observables, from which one can obtain snfficiently inde¬ 
pendent probes simultaneonsly for the compnted initial 
states, for the QCD matter rj/s{T) and also for the appli¬ 
cability of the fluid-dynamical framework: especially, the 
measured centrality systematics of the probability dis¬ 
tributions of 5v2 test the computed initial states, while 
the LHC and RHIC flow-coefficient systematics together 
with the LHC event-plane angle correlations constrain 
the rj/s{T) remarkably consistently. 


II. FLUID DYNAMICS 


Fluid dynamics emerges as an approximation to the 
spacetime evolntion of the system when the microscopic 
scales are small compared to the macroscopic scales like 
the size of the system. Basic equations for flnid dynam¬ 
ics are the conservation laws = 0, and = 

0, where is the energy-momentnm tensor and Nf 
are the possible additional conserved currents (charge, 
baryon nnmber, particle nnmber, etc). In general, 
and can be decomposed w.r.t. the fluid 4-velocity u^, 
defined in the Landau frame eu^ = T^’^Uy, as 

T^^'' = eu^^u'' - ( 1 ) 

Nf = niU^ + n^, ( 2 ) 

where e = T^'^u^Ui, is the local energy density, P = 
Po + n is the isotropic pressnre (snm of eqnilibrinm pres¬ 
sure Pq and bnlk viscons pressnre H), is 

the shear-stress tensor, rii = Nf are the local parti¬ 
cle densities, and nf = are the particle diffusion 

currents. The angnlar brackets indicate the projection 
operators that take the symmetric and traceless part of 
the tensor that is orthogonal to the fluid velocity, i.e., 
and 
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+ A^A); - -A^‘'A„^ 


(3) 


where A^*^ = g^'' — u^u'', and is the metric tensor 
for which we nse the = diag(-|-, —, —, —) convention. 
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The conservation laws are completely general. How¬ 
ever, they are not enough to solve the evolution of the sys¬ 
tem, but additional constraints are needed. In the fluid 
dynamical approximation these additional constraints 
are provided by the evolution equations for the dissi¬ 
pative quantities like For example, in the Navier- 

Stokes (NS) approximation the dissipative quantities are 
directly proportional to the gradients of the equilib¬ 
rium fields (like temperature T, and fluid velocity), e.g., 
= 277(r,{^J)V<'^u'^> and Hns = -C(T, {/r,})V^u^ 
where The microscopic properties of the 

matter are then integrated into the coefficients rj{T, {fJ-i}) 
and Ci'Tj {^i})) which in general depend on the tempera¬ 
ture T and the chemical potentials {^i} associated with 
the conserved charges. It is, however, known that the 
relativistic NS theory is not intrinsically stable, i.e., even 
the hydrostatic equilibrium is linearly unstable |fl5l Ifl6| . 
Therefore, the relativistic NS theory is not suitable for 
the full dynamical description of the system. 


A. Transient fluid dynamics 

The reason for the instability of the NS theory can be 
traced to the fact that the resulting equations of motion 
are parabolic. Therefore, in this theory the signal propa¬ 
gation speed is not limited, and can exceed the speed of 
light, rendering the theory acausal, which in turn makes 
the theory unstable |66|- This problem is solved in the 
Israel-Stewart theory m by taking into account a part 
of the microscopic transient dynamics, e.g. the shear- 
stress tensor relaxes towards the NS values within the 
relaxation time and not instantaneously like in the NS 
theory. The relaxation times are fundamental proper¬ 
ties of the matter similarly to the transport coefficients 
introduced above, and in general they can depend on 
temperature and chemical potentials. 

In this work, we use the equations of motion (e.o.m.) 
derived from kinetic theory EZHZa. Transient fluid dy¬ 
namics can be derived from a microscopic theory by ex¬ 
panding around an equilibrium state and neglecting all 
the microscopic time scales except the slowest one iza. 
This procedure leads to relaxation type equations of mo¬ 
tion for the dissipative quantities, e.g. the evolution equa¬ 
tions for the shear-stress tensor read Eaizi], 

dr 

where the terms up to the first order in gradients (or 
Knudsen number, a ratio of microscopic and macroscopic 
time/length scales, such as Kn ~ iza), second 

order in inverse Reynolds number ^ and product 

of inverse Reynolds and Knudsen number are included. 
Here and i is the 

vorticity tensor. For the purposes of this work, we shall 
neglect the effects of bulk viscous pressure and diffusion 


currents, i.e., H = 0 = nf^. Thus, all dissipative ef¬ 
fects originate in this work only from the dynamics of the 
shear-stress tensor. If one includes also the bulk viscosity, 
several new terms that couple the shear-stress tensor and 
bulk pressure appear also in the e.o.m. of the shear-stress 
tensor |69l [76] . The bulk viscosity can still be important 
around the phase-transition, even if the bulk viscosity is 
negligible in the QGP and the low-temperature hadronic 
phase. However, the magnitude and importance of a pos¬ 
sible large bulk viscosity near the QCD phase transition 
has not yet been fully established [5^ [77H55] . 

Besides affecting the spacetime evolution of the densi¬ 
ties and velocity, viscosity also modifies the local particle 
distributions. For example, in the original work by Israel 
and Stewart m transient fluid dynamics was derived 
from the Boltzmann equation by using the so-called 14- 
moment approximation, where the distribution function 
due to the non-zero shear-stress tensor is written as 


fi(x,p) = foi{x,p) + 6fi = foi[x,p) 


' PrtiPiuT^^'' 

^2T2(e + Po) 


(5) 

Here pf is the 4-momentum of the particle and /oi is the 
equilibrium distribution function. 


foi {x,p) = 


( 2 ^)= 


exp 


T 


Pi 


± 1 


-1 


( 6 ) 


where gi is the degeneracy factor of hadron i. This form 
of Sf does not follow uniquely from the Boltzmann equa¬ 
tion, but is rather the first term of the full moment ex¬ 
pansion E^. Nevertheless, most studies of relativistic 
heavy-ion collisions use this form, and also we adopt this 
procedure here. Currently, the momentum dependence of 
the Sf corrections remains one of the major uncertain¬ 
ties in the fluid dynamical models, see e.g. Refs. |84I[55] 
for studies of the effects of different forms of Sf. For an 
approach to derive Sf corrections from a simplified mi¬ 
croscopic theory, i.e., relaxation time approximation to 
the Boltzmann equation, see EZIlHEj- 


B. Applicability of fluid dynamics 

Fluid dynamics becomes a good approximation when 
gradients are sufficiently small and the evolution of the 
macroscopic variables is slow compared to the micro¬ 
scopic time scales. The systems formed in heavy-ion col¬ 
lisions are, however, very small and their lifetime is short, 
and these conditions are not trivially fulfilled. The esti¬ 
mates of the Knudsen numbers, i.e. ratio of microscopic 
and macroscopic scales, reached in the collisions indicate 
that even with small values of shear viscosity, there can 
still be large corrections to the fluid dynamical evolu¬ 
tion m- Especially in the low density hadronic matter, 
where viscosity is expected to become large EZlIMia, 
the fluid dynamical treatment becomes less reliable. In 
particular, this is true for the decoupling from a fluid to 
free particles, a process that cannot even in principle be 
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fully described by fluid dynamics. Therefore, even if the 
fluid dynamical models have been very successful in de¬ 
scribing the low-pT hadron spectra measured at RHIC 
and LHC energies, it is still not clear in how detail one 
should trust the fluid dynamical description, and what 
are its limitations. 

It is then clear that reaching the final goal of deter¬ 
mining the transport properties of the matter from the 
experimental data requires that also the uncertainties re¬ 
lated to the fluid dynamical evolution are systematically 
charted. There are currently a few ways of extending the 
applicability of fluid dynamics. For example, the moment 
expansion of the Boltzmann equation provides a way to 
include in principle arbitrary orders of the gradients into 
the description, and it has been shown that including all 
the second order terms consistently into the description 
is essential in describing the detailed structure of shock 
waves [93]. One of the characteristics of heavy-ion col¬ 
lisions is that the early expansion is highly asymmetric, 
i.e. the system starts with a fast longitudinal expansion, 
and transverse expansion develops only later. This kind 
of anisotropic expansion results in also highly anisotropic 
local momentum distributions, which can lead to a break¬ 
ing of the usual fluid dynamical description. This is 
the motivation for the so-called anisotropic hydrodynam¬ 
ics |94fl55] . where the functional form of the expansion 
around the equilibrium state is designed to allow large de¬ 
viations from an isotropic momentum distributions. Nei¬ 
ther of these methods are, however, applied to a full de¬ 
scription of heavy-ion collisions, yet. 

One of the important conditions for the applicability 
of fluid dynamics is that different systems should be de¬ 
scribed by the same transport coefficients that can de¬ 
pend on temperature and chemical potentials, but not 
e.g. on the collision energy or the nuclear mass number. 


C. Our fluid dynamical setup 


In this work we employ the setup previously used 
in Refs. m m m [^, where the longitudinal ex¬ 
pansion is approximated by a scaling flow consistent 
with longitudinal boost-invariance. In this approxima¬ 
tion the longitudinal flow velocity is given by = z/t, 
and the components of the energy-momentum tensor, 
Eq. Q, become independent of the spacetime rapidity 
r]s = ( 1 / 2 ) In [(t-I- 2 :)/(t — z)], i.e., they depend on the 
transverse coordinates, r = {x,y), and the longitudinal 
proper time, r = y/t'^ — z^, only. From a numerical point 
of view, this reduces the (3+l)-dimensional problem to 
a ( 2 +l)-dimensional one. 

The coefficients of the non-linear terms in the equa¬ 
tions of motion for the shear-stress tensor, Eq. ([4]), are 
taken from the 14-moment approximation to the ultra- 
relativistic gas jSS] ISni Elj, i.e.. Cl = — (4/3)t^, C 2 = 
— (10/7)r^, C3 = 2r^, and C4 = 9/(70Po), and the relation 



FIG. 1. (Color online) Parametrizations of the temperature 
dependence of the shear-viscosity to entropy ratio, labelled 
here in the order of increasing rj/s at T = 100 MeV. For more 
details, see the text and Table |T| 


between the relaxation time T,r and the shear viscosity is 


Tn 


5r] 

e + Pq 


(7) 


In thermodynamical equilibrium, the properties of the 
matter are essentially given by the EoS that gives pres¬ 
sure as a function of temperature. Here we use the 
s95p-PCE-vl parametrization of lattice QCD results at 
zero net-baryon density EZj. The high-temperature part 
of this EoS is from the hotQCD collaboration |98l 
and it is smoothly connected to a hadron resonance gas, 
where resonances up to mass of 2 GeV are included. The 
hadronic part of the EoS includes a chemical freeze-out 
at Tchem = 175 MeV, where all stable hadron ratios are 
fixed 110014102] . A hadron is considered stable, if its life¬ 
time is more than 10 fm. In the perfect fluid limit the 
construction of the chemical freeze-out also conserves the 
number of stable particles. However, in the viscous fluid 
there is still small (approximately 1 %) entropy produc¬ 
tion below Tchem = 175 MeV, and this leads to a small 
increase in the number of particles during the evolution 
of chemically frozen hadronic matter. 

Once the transport coefficients and EoS above are 
given, the only degrees of freedom left are the shear vis¬ 
cosity to entropy density ratio rj/s{T) and the initial com¬ 
ponents T'^'^(ro, r). In the boost-invariant approximation 
it is enough to specify r) in the transverse plane 

at some initial proper time tq. The initial conditions 
calculated from the EbyE EKRT setup are discussed in 
detail in the next section. 

As shown in Fig. we parametrize the temperature 
dependence of the r]/s ratio in a similar manner as we did 
in 1 5 5] . by assuming a minimum of 77/5 at T = T„iin to 
be somewhere in the cross-over temperature-region and a 
linearly rising (decreasing) behavior in the QGP (HRG) 
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phase. Table l^shows the corresponding parameters from 
which these linear slopes can be constructed. We have 
converged into these parametrizations iteratively, requir¬ 
ing them to reproduce the measured 2-particle cumulant 
elliptic flow V2{2} (see Sec. IVB for the definition) in 
mid-peripheral collisions at the LHC. In addition, we also 
exploit the HH-HQ parametrization of Ref. [Hiiiiim] 
(used later also in Ref. |7]), which features a rapid growth 
of r]/s(T) in the QGP combined with a more modest 
decrease in the hadron gas phase. We label the above 
parametrizations here as paraml, param2, paramZ = 
HH-HQ, and paramA, in the order of an increasing value 
of 77 /s at Tdec = 100 MeV. As we will show, a simul¬ 
taneous comparison with the RHIC results is then nec¬ 
essary to see the sensitivity to ri/s{T). As indicated in 
Fig. [3 we perform the calculations also for a constant 
rj/s = 0.2, keeping also this value unchanged from the 
LHC to RHIC. The sensitivity of the computed Vn to a 
constant rj/s = 0.2 ±0.1 will be demonstrated. 


TABLE I. The constant-slope parametrizations of rj/s{T), 
constructed so that they reproduce the LHC Vn data. 

Tmin/MeV {rj/s)m\n ?7/s(100 MeV) 77/s(500MeV) 


paraml 

150 

0.12 

0.24 

0.65 

param2 

180 

0.16 

0.36 

0.16 

param4 

180 

0.12 

0.76 

0.30 


Once the initial conditions, EoS, and the transport co¬ 
efficients are given, the equations of motion for shear- 
stress tensor, Eq. Q, and the conservation laws form a 
closed system of equations that can be solved numerically 
to obtain the spacetime evolution of all the quantities 
appearing in the energy-momentum tensor, Eq. Q . The 
numerical algorithm employed here to solve the equations 
of motion is introduced and discussed in Refs. mm- 


In this work we take the freeze-out surface to be 
a constant-temperature surface with Tdec = 100 MeV, 
which gives a good agreement with the slopes of the mea¬ 
sured charged hadron pT spectra. A more physical way 
would be to decouple the system dynamically on a sur¬ 
face where the expansion rate of the system becomes of 
the same magnitude as the average scattering or thermal- 
ization rate (here T,r), he. when Kn ~ 1 |251 llOflHlll] . 
However, in practice, the differences to the constant- 
temperature freeze-out are quite modest, especially near 
midrapidity. 

In principle, the Cooper-Frye integral ([8 1 should be 
calculated for all the hadronic states included into the 
EoS, i.e., up to a mass 2 GeV. However, in order to save 
computational time, we include here hadrons only up to 
a mass 1.5 GeV. In practice, the effect on the final results 
shown here is negligible. All the strong and electromag¬ 
netic two- and three-particle decays of the hadronic res¬ 
onances (most of the hadrons in the EoS are unstable and 
decay before they can be observed) are calculated here ac¬ 
cording to Ref. In finding the constant temperature 

hypersurfaces, we employ the Cornelius algorithm |113| . 


III. INITIAL CONDITIONS FROM THE LOCAL 
EKRT SATURATION MODEL 

Let us then discuss the details of the NLO-improved 
pQCD ± local saturation framework [551155] , which com¬ 
bines a NLO pQCD computation of the minijet trans¬ 
verse energy Et production with saturation of gluon pro¬ 
duction. First, we discuss the computation for averaged 
(smooth) initial conditions, after which we explain how 
the event-by-event setup utilizes these calculations. 


A. Minijet Et production in A±A collisions 


D. The freeze-out stage 

The fluid dynamical quantities are not directly com¬ 
parable to the experimental data. Therefore, it is nec¬ 
essary to convert them into the experimentally observ¬ 
able hadron transverse momentum spectra. Here we em¬ 
ploy the standard Cooper-Frye procedure una, where 
the spectrum is calculated as the number of particles 
crossing some surface S whose normal vector is 
This leads to a Lorentz-invariant spectrum for a hadron 


_ 1 

d^p 2 dydp\,d(j) 


d^^f,ix)p^f{x,p), ( 8 ) 


where p^ = {E, p) denotes the four-momentum of the 
hadron, and f{x,p) is the single-particle distribution 
function, Eq. (j^), of the hadron on the surface. In the 
boost-invariant approximation the spectrum is indepen¬ 
dent of the rapidity y. 


For a given collision energy ■sJsnn and nuclear mass 
number A the initial minijet Ej- produced perturbatively 
into a rapidity window Ay in A+A collisions and above 
a transverse momentum scale po ^ Aqcdj can be com¬ 
puted as [53] 

d-E'j’ 

y/sNN, A, r, b; jd) = Ta{vi)Ta{v2)(t{Et,) poAyA 

(9) 

where ri/2 = r ± b /2 with r = {x, y) denoting the trans¬ 
verse coordinate and b the impact parameter. The nu¬ 
clear collision geometry is given by the nuclear thickness 
function 

/ OO 

dzpA{r,z), (10) 

-OO 


where the nuclear density pa{x,z) is parametrized with 
the standard Woods-Saxon (WS) profile 


PA{r,z) 


no 


exp 




±1 


( 11 ) 
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with the nuclear radius Ra = —0.86A“^/^) fm, 

d = 0.54 fm, and ng = 3A/(47ri?^)[(l + T:'^d?/R\)]~^ « 
0.17 fm“^. According to collinear factorization and 
pQCD, the first A-r-nioment of the minijet Et distribu¬ 
tion, (j{ET)po,Ay,i 3 , in NLO is computed as [571 IITI155] 

(^{Et)p„,Av,p = / AEtEt-^ , ( 12 ) 

Jo (IEt pg,Ay,p 

where the semi-inclusive Et distribution of minijets in a 
rapidity interval Ay in N+N collisions is given by 


dcr 


dEj 


^ 1 


PoAy,P n=2 


dcr^-*'" 


(13) 


Here, the n-particle momentum phase-space integration 
d[PS].,j takes place in 4 — 2£ dimensions, and we have 
introduced a compact notation for the differential NLO 
partonic cross sections dcr^“*'^/d[PS ]2 and dcr^^^^/dpSJg, 
corresponding to the (2 —)■ 2) and (2 —)■ 3) scatterings, 
respectively. A detailed discussion of the d(T^“*'”/d[PS]^ 
which consist of (NLO, MS scheme) PDFs and squared 
spin- and color-summed/averaged scattering matrix ele¬ 
ments, summed over all possible parton types, is given in 

mm- 

The IR and CL singularities present in the partonic 
cross sections at order are regulated by computing the 
(2 —)■ 2) and (2 —)■ 3) squared matrix elements in 4 — 2e 
dimensions. The ultraviolet divergences present in the 
(2 —>■ 2) parts are taken care of by renormalization using 
dimensional regularization and the MS scheme. The full 
analytical calculation for these squared matrix elements 
was done first in |115| . and details of some of these rather 
complicated calculations are given in [116] . The phase 
space differentials d[PS ]2 and dpSJg stand for 

d[PS]2 = dpT2dj/idj/2d^"^"())2, 

dpSJg = dpT2dpT3d2/idy2d2/3d^“^"(/>2d^“^''<?i3, 


where the appropriate kinematical variables for the two- 
and three-parton phase spaces are the transverse mo¬ 
menta pTi = |pTi|j rapidities yt and azimuth angles 
(f>i. For the two-parton final state, the transverse mo¬ 
mentum conservation determines pxi = pt 2 and (j)i = 
(j )2 + TTj and similarly for the three-parton final state 
Pti = — (pr 2 + Pts)- The measurement functions S 2 
and iSg in Eq. (131 specify the physical quantity to be 
computed. As explained in na, the cancellation of the 
remaining IR and CL singularities between the UV renor¬ 
malized squared (2 —>■ 2) and (2 —>■ 3) matrix elements 
takes place only if the measurement function S'g reduces 
to the S 2 in the soft (the energy of one of the final-state 
partons vanishes) and collinear (one of the final state par¬ 
ticles becomes collinear with any other particle) limits. 

In our case, the measurement functions define the to¬ 
tal minijet Et produced into a mid-rapidity window Ay 
defined in the (y, ^)-plane as 


The minijet Et entering Ay is defined here as a sum of 
the transverse momenta pTi of those final-state partons 
whose rapidities are in Ay 

n—2,3 

E'j' = E Q{y^ e Ay)pTi, (16) 

i=l 

where all partons are assumed massless and 0 is the 
standard step function. For computing the minijet Et 
distribution, our measurement functions must also spec¬ 
ify which scatterings are to be considered hard and 
thus included in the perturbative calculation. We de¬ 
fine the hard perturbative scatterings to be those with 
large enough transverse momentum produced, regardless 
of where the partons go in rapidity, 

n=2,3 

E PTi > 2po, (17) 

1=1 


where po > Aqcd- 

Now, for the (2 —>■ 2) hard processes transverse mo¬ 
mentum conservation ensures that if at least one parton 
falls into our rapidity acceptance, then Et > Po- How¬ 
ever, in the (2 —> 3) case we may have processes which ful¬ 
fil the requirement of being hard [pti +Pt 2 +Pt 3 > 2po) 
but bring less than po of Et in Ay. This happens, e.g., 
for configurations where two hard partons fall outside 
Ay and only one softer parton with pT < po enters Ay. 
Therefore, the remaining freedom in defining the mea¬ 
surement function S 3 is that in the (2 —> 3) case we may 
still restrict the amount of the minimum Et in Ay in 
an IR/CL safe way. In [53] it was shown that in the 
S 3 case in fact any minimum amount, Et > Ppo, where 
0 < ,5 < 1, gives an equally good IR/CL safe restriction 
for the Et in Ay, which relaxes back to the S 2 case in 
the soft and collinear limits. 

Thus, the IR- and CL-safe measurement functions S 2 
and S 3 can now be written down by combining the def¬ 
inition of minijet Et in Ay, the definition of the hard 
perturbative scatterings and the restriction of minimum 
Et discussed above. 


= (5 Et — 


'^e{yiG Ay)pTi 




(18) 


X 0 ( Et"^* - 1 ^ > Ppo), 


\i=l 


where /3 is a phenomenological parameter to be deter¬ 
mined from the experimental data. Next, integrating the 
delta functions away in Eq. (121 we obtain 


da^ 


0’(£’r)po,Ay,/3 — E/ yci[PS] 


• I 2 /I < 0.5, 0 < (/ < 27r. 


(15) where the IR and CL safe measurement functions for the 




















first Et moment are denoted by 


Sn = Q{y^ e ^y)PT^ j x 0 PTi > 2po 



x0 


^0(?/j e Ay)pTi 


( 20 ) 


>(3po]- 


The numerical computation for the rather complicated 
six-dimensional integrals uni in Eq. ( |19[ ) is performed 
with Monte Carlo integration, using an updated version 
of the code developed for Uni im 1^ where the (2 —)• 3) 
parts and their partonic book-keeping are based on the 
Ellis-Kunszt-Soper jet code |118| |114j . For the DGLAP 
evolved nPDFs, we apply the NLO CTEQ6M free pro¬ 
ton PDFs |119| together with the latest set of transverse- 
coordinate (Tyi(r)) dependent NLO EPS09s nuclear ef¬ 
fects ino]. The implementation of these spatial nuclear 
effects is done as instructed in m, calculating the re¬ 
sults directly for each r and b; for details, see |12()| . The 
renormalization scale pn and factorization scale pp are 
chosen equal, pR = pp = P- We set the scale p to be 
proportional to the total transverse momentum produced 
in the hard perturbative scattering, regardless of the par- 
tons being in Ay or not: 


C 


= v S 


PTi 


( 21 ) 


where the constant C is set to unity. 


B. Local saturation of minijet Ep production 

As explained in |53| . the low-transverse-momentum 
parton (dominantly gluon) production can be conjec¬ 
tured to be controlled by saturation of minijet Ep pro¬ 
duction. In this new EKRT approach the saturation 
takes place when (3 —>■ 2) and higher-order partonic pro¬ 
cesses start to dominate over the conventional (2 2) 

processes (and (2 3) at higher orders). Thus, at satu¬ 

ration, we require that the rapidity densities of the pro¬ 
duced Ep fulfil the condition 




( 22 ) 


To LO in as, the l.h.s. scales as 


ji^(2^2)~(r«)^(|)„, (23) 

where we assign the factor Tap for each of the incoming 
gluons, a^lpQ for the a{2 —>■ 2) partonic cross section and 
the cut-off scale po for the Ep. Here, g denotes the gluon 
PDFs. Similarly, for the r.h.s (3 —2) term in Eq. (22) 
we may write 

1^(3 ^ 2) ^ (TAgf\ f 4) PO, (24) 

d^rdy pg \pf^ J 


where the scale Pq ^ is to compensate the fm ^ dimension 
of the extra Ta in Eq. (24). Substituting the Eqs. (23) 


and (24) into the saturation condition (22), we get 


(r.9)=(f)p.~(Wi(f)p., (25) 

which leads to a scaling Tap ^ Polctg for the gluon 
density probed at saturation. Feeding this scaling law 
back to the saturation condition in Eq. (22), we obtain 


a transversally local saturation criterion for the minijet 
Ep production in A+A collisions at non-zero impact pa¬ 
rameters 1551. 


d-Ep 

d^r 


(po, y/sNN, A, Ay, r, b; /3) = ^^p^Ay, (26) 


with an unknown (but to a first approximation ag- 
independent) proportionality constant ATgat 1, whose 
value needs to be determined from the data. Once 
the saturation scale is obtained as the solution po = 
Psat(V^JVAT, A, Ay, r, b; /3, ATg^t) of Eq. ([^, we get the 
total amount of minijet transverse energy dEp{pQ = 
Psat)/d^r produced into a mid-rapidity window Ay. 


C. Numerical implementation 

The procedure to obtain the locally saturated NLO 
minijet Ep is straightforward, but the challenges in the 
numerical implementation are worth mentioning. First, 
with the spatially dependent nPDFs the computation of 
the locally saturated NLO dAT/d^r(psat) becomes slow, 
mainly due to the multidimensional MC integrations in 
the (2 —>■ 3) parts. Second, since pgat can be deter¬ 
mined from Eq. (26) through iteration only, we need 
dAj’/d^r(po) for 0(10) different po’s at each r for each 
b. Third, the spatial {x,y) grid for constructing initial 
conditions for fluid dynamics has to be dense enough, 
say Ax = Ay = 0.4 fm, and extend far enough, at least 
to r ~ Ra where the approach can still be imagined to 
work. In one quarter-plane we then have to compute the 
saturated minijet Ep in 0(250) different {x,y) points for 
each b. Fourth, and worst, we have to determine the free 
parameters ATgat and /3 iteratively on the basis of the cen¬ 
trality dependence of the bulk data, i.e., after performing 
the hydrodynamic evolution for all centrality classes with 
initial conditions computed for each ATgat, j3 pair with a 
given g/s. Thus, a blindly repeated NLO computation 
of locally saturated averaged initial conditions for such 
an iterative procedure becomes numerically too slow, and 
the EbyE framework would then seem just impossible. 

The first key-observation in circumventing the above 
critical slowness problems, made in [^, is that to a good 
approximation the "K-factor" 


AT = cr(AT)po,Ay,/3(NLO)/cr(£’T)po,Ay,/3 (lo) (27) 

does not depend on the PDFs (free proton, nuclear or 
spatial). Then, the full NLO result can be approximated 
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FIG. 2. (Color online) Saturation momentum paat as a func¬ 
tion of nuclear overlap density TaTa with ifsat = 0.5 (a) 
and Kaat = 0.75 (b) in Pb+Pb collisions at the LHC (red 
points), and in Au+Au collisions at RHIC (blue points), cal¬ 
culated with several different impact parameters. The dashed 
lines show the corresponding parametrization | |29[ ), and its ex¬ 
trapolation to the typical highest TaTa’s we encounter in the 
EbyE analysis. 


nuclear overlap density, 


Paa (r) = (^r - Ta , (28) 

in y/sNN = 200 GeV Au+Au collisions at RHIC and in 
= 2.76 TeV Pb+Pb collisions at the LHC with 
P = 0.8 and ATgat = 0.5 (a) and 0.75 (b). The blue and 
red points in the figure are from the pQCD+saturation 
calculation at different transverse positions and with sev¬ 
eral different impact parameters. Thus, for a fixed cms- 
energy and collision system, the computed local satura¬ 
tion scale Psa,t{x, y) is to a very good approximation only 
a function of the paa, and furthermore, the function is 
the same for all centrality classes. 

The emergence of such a scaling can be understood 
as follows. In the naive scaling limit, where the mini¬ 
jet a{ET) oc Pq^, the saturation criterion (261 leads to 
the scaling p'^^^ oc {paaY with 5 — 1/2. Ai'^discussed 
in Ref. [M] (in LO, without nPDFs), corrections to the 
power 5 can be traced back to the x— and Q^-slopes of 
the small-a; gluon distribution, phase-space integration 
and running of a^. 

Figurej^now shows that also the NLO calculation with 
nPDFs preserves the power-law scaling property of Psat 
extremely well for a fixed cms-energy and for a fixed nu¬ 
cleus A. The spatial effects in the nPDFs could still mod¬ 
ify this scaling from one impact parameter to another. 
Figure shows, however, that the these effects are so 
small that the paa dependence of p^^t is to a good ap¬ 
proximation universal over all centralities. Thus, we can 
very accurately parametrize the saturation scale as 


PsatipAA) = C [a + paaT - bCa^ 


(29) 


where a, b, C and n are parameters that depend on 
A, y/sNN, ATsat and /3. For a given A and -^/snn the 
(A^sat,/3)~dependence can be parametrized by a polyno¬ 
mial. 


T)(A"sat, P) — + anKsat + ^i2P 

+ aaKsatP + anP"^ + 


(30) 


by implementing the spatial nPDFs into the fast LO part 
only, and using the K-factors to account for the NLO 
effects, i.e. 


Cr(£''r)po,Ay,/3(NLO,EPS09s) « cr(£’T)po, Ay, /3(lO,EPS 09s) X AT, 


The coefficients Oy for the parameters a, b, C and n 
are listed in Tables |H| |HI| I^ and|V|for = 2.76 

TeV Pb+Pb collisions and y/sNN = 200 GeV Au+Au 
collisions. Note that the parametrizations are found sep¬ 
arately for P < 0.9 and P > 0.9. Armed with the above 
parametrization of PsatiPAA), we have been able to chart 
the Ksat,P plane for finding the initial conditions dis¬ 
cussed next, and develop the EbyE framework. 


where the K-factor has been computed only once, with 
the free-proton PDFs. According to the checks we have 
made over the (x, y) plane, this approximates the full 
NLO result very well, within a few percents both at RHIC 
and LHC. 

The second key-observation enabling the locally satu¬ 
rated EKRT framework is demonstrated in Fig. which 
shows the calculated values of psat as a function of the 


D. Initial state for fluid dynamical evolution 

As initial conditions, our boost-invariant dissipative 
fluid-dynamical modeling requires the transverse energy 
density e(r, tq), transverse velocity VT(r,ro) and initial 
shear stress tensor 7r^'^(r,ro) at a constant initialization 
proper time tq of fluid-dynamics. 
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TABLE II. The parametrization of Psat(A'sat, /3) for ^/snn = TABLE IV. The parametrization of Psat{Ksa.t, /3) for y'siviv = 
2.76 TeV Pb+Pb collisions for iLsat £ [0.4, 2.0] and j3 < 0.9 200 GeV Au+Au collisions for iLaat £ [0.4, 2.0] and /3 < 0.9 


C 


Oio 3.9027590 0.1312476 
an -0.6277216 -0.0157637 
aa 1.0703962 -0.0362980 
flis 0.0692793 -0.0022506 
an -1.9808449 0.0615129 
Ois 0.1106879 0.0052116 


-0.0044020 0.8537670 
0.0220154 -0.0580163 
-0.0005974 0.0957157 
0.0125320 -0.0016413 
-0.0032844 -0.1788390 
-0.0033841 0.0220187 


P^ 


C 


aio 

an 

ai2 

ai3 

an 

ais 


10.3313939 0.0303079 -0.0070317 0.9381026 
-0.3165983 -0.0024562 0.1561924 -0.0005718 
-12.8128174 0.0139955 -0.0026174 0.0376918 
-0.0273664 -0.0017971 -0.0369552 0.0072667 
4.6810067 0.0923750 -0.0174187 -0.3018326 
0.0527041 0.0005875 -0.0226980 0.0013976 


TABLE III. The parametrization of psat (A'sat, /3) for -sJsnn = TABLE V. The parametrization of psat (Tfsat, /?) for ^/snn = 
2.76 TeV Pb+Pb coiiisions for Aiaat £ [0.4, 2.0] and /3 > 0.9 200 GeV An+An coiiisions for iLaat £ [0.4, 2.0] and P > 0.9 


Pi 

aio 

an 

an 

an 

an 

an 


C 


n 


a 


27.3259359 -1.9924684 0.1038047 0.5211725 
-0.3371381 0.0835716 0.0539039 -0.6286044 
-42.6176287 4.1698751 -0.2099840 2.5059182 
-0.1844621 -0.1206132 -0.0144174 0.7131778 
17.6786774 -1.9891770 0.0950212 -2.5125962 
0.3092463 0.0003279 0.0014117 0.0150475 


C 


n 


a 


ffliO 

an 

an 

an 

ai4 

an 


91.4314177 -0.4406026 0.7332375 3.0875818 
2.5123667 0.0782859 0.2132747 -0.2205018 

-165.8206094 0.6486681 -1.5009886 -3.8563125 
-2.6487281 -0.1005554 -0.0219393 0.2777689 
77.0170469 -0.0909378 0.7419402 1.5327054 
0.2192064 0.0004503 -0.0336409 -0.0006138 


In this work the initial transverse velocity and shear 
stress tensor are chosen to be zero. The transverse pro¬ 
file for the local initial energy density at the formation 
(production) of the system is computed similarly as in 
Refs. |5^HE1I55] . 


e(r,Ta(r)) 


d^T 1 
d^r Ts(r)Ay 


-Rsat 

TT 


[Psat(r)]'‘, 


(31) 


where the local formation time of the minijet plasma at 
each transverse point r is given by Ts(r) = l/psat(r). 
Since for the fluid-dynamical evolution we need the initial 
state at a fixed time, the computed energy densities have 
to be evolved to the same tq at each r. To do this, we first 
set a minimum scale = 1 GeV for which we assume 
that we can still trust the pQCD calculation. This corre¬ 
sponds to a maximum formation time tq = l/p““ « 0.2 
fm in our pQCD+saturation setup. Next, the uncertain¬ 
ties in the "pre-thermal" evolution from Ts(r) to tq can 
be studied by considering the two limits: 1) the Bjorken 
free streaming (FS) scaling 


our final results will be relatively insensitive to the pre- 
thermal evolution. For this reason, in the present study 
we stick to the latter (BJ) case. 

Finally, we need the initial energy densities at the 
edges of the system which are outside the applicabil¬ 
ity region of our pQCD+saturation model, i.e. the en¬ 
ergy densities below Cmm = -K^sat[p™]^ at tq. To obtain 
these, we smoothly connect the BJ-evolved energy den¬ 
sity to the binary profile, i.e. the energy density profile 
is parametrized below Cmin as e = C(TaTa)^, where the 
power n is given by 


1 

^=2 


(fc + 1) + (/c — 1) tanh 


/ ctnnTaTa — g 
\ ^ 


, (34) 


with the total inelastic nucleon-nucleon cross-section unn 
and g = S = 0.5 fm“^. The parameters C and k are 
constants that ensure a smooth connection at e = emin- 


E. Averaged initial conditions 


e(r. To) = e(r, rs(r)) (32) 

which preserves the transverse energy, and 2) the Bjorken 
hydrodynamic scaling solution (BJ) 

e(r,To) = e(r,rs(r)) , (33) 

where a maximum amount of energy is transfered into the 
longitudinal direction by the P^dV work. As discussed in 
|55j . due to the freedom we still have in fixing (iLgat,/^), 


As an example we show the calculated initial energy 
density profiles in y/sNN = 2.76 TeV Pb+Pb collisions 
at To = 0.20 fm in 0 — 5 % and 20 — 30 % centrality 
classes in Figs. and respectively. The calculation 
of the nuclear overlap geometry and of the impact pa¬ 
rameters corresponding to the centrality classes are in 
this case based on the optical Glauber model. For com¬ 
parison, we also show the usual simple Glauber model 
based eBG and eWN profiles |53]. The eBG and eWN 
profiles are normalized such that the initial entropy per 
unit spacetime rapidity, dSi/drjs, which in the ideal fluid 
is directly proportional to the final hadron multiplicity, is 
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the same as in the calculated initial state in 0 — 5 % cen¬ 
trality class. Overall, the energy density gradients from 
the EKRT model are slightly steeper than in the eWN 
profile, but not as steep as in the eBC profile. 

The initial profiles can be further quantified by calcu¬ 
lating the eccentricity, 

= -{r™e“^}/{r™}, (35) 

where the curly brackets denote the average over the 
transverse plane, i.e., {•••} = / dxdy e[x,y,To){- ■ ■), r is 
the distance to the system’s center of mass, and e(x, y, tq) 
is the energy density at the initial time rg. The “partici¬ 
pant plane”-angle 4'm,n can be calculated as 

'km n = —atan2 ({r™ cos(n^)}, {r™ sin(n(/>)}) -I- —, (36) 
’ n n 

where the atan2(x, y) function gives the angle in the 
correct quadrant of the transverse plane. In the ab¬ 
sence of event-by-event fluctuations the event-plane angle 
^m,n = 0, if the a;-axis is chosen in the direction of the 
impact parameter, and em,n = 0 for all odd n. Note, 
however, that later when we consider the event-by-event 
density fluctuations the phase and em,n for odd n are not 
generally zero, but fluctuate from event to event. We also 
use a short-hand notation = £n,n- The eccentricities 
£2 of the calculated initial profiles as a function of cen¬ 
trality are shown in Fig. As before, we show the 
comparison to the eBC and eWN profiles, and we can 
immediately see that £2 of the pQCD-based initial con¬ 
ditions are between the eBC and eWN Glauber model 
limits. 

The corresponding energy density profiles in ■sJsnn = 
200 GeV Au+Au collisions are shown in Figs. and 
and the initial eccentricities in Fig. Overall the 
pQCD initial states are quite similar at RHIC and the 
LHC. The most notable change is that the eWN initial 
state is closer to the pQCD initial state at RHIC energy, 
as can be seen both in the energy density profiles and 
eccentricities. 

The computed energy density profiles discussed above 
were used as initial conditions to fluid dynamical evo¬ 
lution in Ref. jSS]. It was shown that this model can 
reproduce the centrality dependence of the multiplicity, 
PT-spectra and elliptic flow coefficients simultaneously at 
the RHIC and LHC energies. However, in order to com¬ 
pare to the available experimental data in more detail, it 
is necessary to take into account the event-by-event na¬ 
ture of the collisions. Inclusion of the effects of the den¬ 
sity fluctuations to the pQCD initial state is described 
next. 


F. Event-by-event density fluctuation 

The main source that drives the initial state density 
fluctuations are the random fluctuations in the positions 
of the nucleons inside the colliding nuclei. Therefore, 


the basic ingredient in modeling such fluctuations is the 
spatial distribution of nucleons inside the nuclei. These 
distributions are mainly constrained by the measured nu¬ 
clear charge distributions. The nuclear charge density 
is frequently parametrized by the Woods-Saxon function 
(111, with R and d as the free parameters. However, 
the measured charge distribution is not the same as the 
nucleon position distribution, because the nucleons are 
not point-like particles, but have a finite size and charge 
distributions themselves. Thus, in principle, the Woods- 
Saxon parametrization for nucleon position should be 
constrained in such way that when folding with the nu¬ 
cleon charge profile it gives the measured nuclear charge 
distribution. The situation is complicated even more by 
the fact that protons and neutrons are not distributed in 
the same way, but especially in heavy nuclei the charge- 
neutral neutrons tend to form the outer layer of nuclei. 
The formation of this “neutron skin” should be taken into 
account when constraining the distributions. In our case, 
however, since we are mainly interested in gluons, whose 
distribution in protons and neutrons are similar, only the 
average nucleon distribution matters. 


Here, in building the EbyE setup, we take the nucleon 
distribution in a Pb nucleus from Ref. m, which is 
already constrained by the charge distribution and avail¬ 
able measurements of the neutron skin thickness. In 
practice, this nucleon density profile can be parametrized 
by the usual Woods-Saxon function, with R = 6.7 fm 
and d = 0.55 fm. It is noteworthy that the neutron 
skin and the finite size of the nucleons tend to affect the 
parametrizations in opposite direction: the final Woods- 
Saxon parameters are actually within the errors of the pa¬ 
rameters given for the measured charge distribution |122j . 
Therefore, effectively we can take the Woods-Saxon pa¬ 
rameters for the charge density and interpret the result¬ 
ing profile as a nucleon position distribution. The the¬ 
oretical models indicate that the neutron skin thickness 
varies only slowly with the nuclear mass number |123| . 
and therefore we can expect that a similar cancellation 
happens also for other heavy nuclei. For Au nuclei we, 
therefore, take the Woods-Saxon parameters from the 
charge distribution and interpret the resulting distribu¬ 
tion as a nuclear position distribution. 


The nucleon positions inside the nuclei are then sam¬ 
pled according to the Woods-Saxon distribution by as¬ 
suming them uncorrelated, i.e. sampling each nucleon 
position independently. In doing this, we keep in mind, 
however, that in principle the nucleon positions are cor¬ 
related, e.g. two nucleons cannot overlap, but in practice 
the effect of the correlations is rather weak m , ex¬ 
cept perhaps in ultracentral collisions |125| . As a result 
we obtain an ensemble of nuclear configurations char¬ 
acterized by the nucleon positions {xi,yi, Zi). By ran¬ 
domly sampling the impact parameters from a distribu¬ 
tion dN/db^ cx const, we then get an ensemble of nuclear 
collisions. 
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FIG. 3. (Color online) Energy density profiles in y/s^N = 2.76 TeV Pb+Pb collisions at tq = 0.20 fm, in the 0 — 5 % (a) and 
in the 20 — 30 % centrality class (b) computed with TsTsat = 0.63 and /? = 0.8 in the BJ prethermal evolution case. The small 
vertical lines show approximately where the matching to the TaTa profile is done, i.e. at psat = 1 GeV. The figures (c) and (d) 
show the same for y/SNN = 200 GeV Au+Au collisions. 



FIG. 4. (Color online) Initial eccentricity as a function of centrality in ^snn = 2.76 TeV Pb+Pb collisions (a) and in 
y/SNN = 200 GeV Au+Au collisions (b) at to = 0.20 fm, computed with Jfsat = 0.63, /3 = 0.8 and BJ prethermal evolution. 
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G. Nuclear and nucleon overlap densities 


In the EKRT minijet framework a nuclear collision is 
regarded as a collision of two gluon clouds rather than a 
collection of individual nucleon-nucleon collisions. The 
leading idea in our EbyE setup is that we first form 
the nuclear overlap density paa locally in r = {x, y) for 
each nuclear collision event, accounting for the nucleon 
configurations in each collision. Then, the local satura¬ 
tion scales _Psat(pAyi(r)) in each event are obtained from 
Eq. (291. The initial energy densities at fixed tq can then 
be computed, EbyE, as described in Sec. |IIID| 

We define the nuclear thickness function Ta in each 
event as a sum of the corresponding nucleon thickness 
functions Tn, 


T. 4 (r)=^r„(|r-r,|), 


(37) 


where the sum is over the nucleon positions in the nucleus 
A and where the Tn have been normalized to one. The 
nuclear overlap density pab (r) in each AtB collision is 
then obtained from Eq. ( [2^ . 

Since the minijet production considered here is dom¬ 
inated by gluonic channels, the r„ above is to be un¬ 
derstood as the gluonic thickness function rather than 
the one obtained from the (better known) charge den¬ 
sities of nucleons. To obtain the gluonic T„ needed 
here, we exploit exclusive electroproduction of J/ij: at 
HERA, 7 * -f p —>■ J/'f/j -I- p, for which ZEUS has mea¬ 
sured the differential cross section near t = 0 to be 
da/dt (X exp(—&|t|) oc \G\^ with a slope b = 4.72 GeV“^ 
|63| . Taking a 2-dimensional Fourier transformation of 
the corresponding 2-gluon form factor G leads to a Gaus¬ 
sian distribution for Tn, 

where the width parameter a = Vb « 0.43 fm. 


H. Centrality selection and sampling the nuclear 
collisions 

After sampling the nucleon configurations and the im¬ 
pact parameter we determine whether a nuclear collision 
occurs by using the following geometric collision crite¬ 
rion: the A+B collision takes place if the transverse dis¬ 
tance between at least one of the nucleons from A and 
one from B is shorter than where unn is the 

total inelastic NN cross-section. At the LHG ctnn = 64 
mb and at RHIG (Tnn = 42 mb. We emphasize that unn 
is here only used in the above collision trigger criterion, 
and that the calculation of the initial state is otherwise 
essentially independent of (Tnn- Following this procedure, 
we create a large number of nuclear collision events, for 
which we then calculate the initial energy density profiles 
as described in the previous sections. 



FIG. 5. (Color online) Probability distribution of the charged 
hadron multiplicity dN^i/drip for the five different y/s{T) 
cases of Fig.[2 compared with the parametrization of the AL¬ 
ICE VZERO amplitude read off from Ref. |126| . in ^snn = 
2.76 TeV Pb+Pb collisions. 


Next, the fluid-dynamical evolution is calculated sep¬ 
arately for each event, after which we calculate the px- 
spectrum and multiplicities as described in Sec. |IID| The 
events are then divided into centrality classes according 
to their final multiplicity (or equivalently the final total 
entropy). For example, the 0 — 5 % centrality class con¬ 
sists of the events with the highest multiplicity, the top 
5 % of the total number of events. 

In Fig. [^we show the calculated probability distribu¬ 
tion of the charged hadron multiplicity dN^h/dp-p com¬ 
pared to the parametrization of the ALIGE measurement 
of VZERO amplitude, read off from Fig. 10 of Ref. |126| . 
that is approximately proportional to the final state mul¬ 
tiplicity. The distributions are scaled to have approxi¬ 
mately the same average. As one can see from the figure, 
the agreement between our calculation and the ALIGE 
measurement is very good, except in the very central 
collisions. This is, indeed, expected as this tail of the 
distribution is dominated by the dynamical multiplicity 
fluctuations which we do not yet include in the current 
EKRT framework. In our case, such dynamical fluctua¬ 
tions would mean that for the same value of the overlap 
density paa = TaTa the saturation scale (i.e. gluon 
multiplicity), and hence entropy production, would be 
fluctuating from one event to another. 

Even without the fluid-dynamical evolution, it is possi¬ 
ble to estimate the centrality dependence of multiplicity 
from the initial entropy. Because the hadron multiplicity 
is proportional to the final entropy (by a factor that de¬ 
pends on the decoupling temperature), the entropy pro¬ 
duction during the fluid dynamical evolution can signif¬ 
icantly affect the multiplicity, but the its effect on the 
relative centrality dependence of the multiplicity is much 
weaker, see below. Once we have the energy density pro¬ 
files, we can convert them to the entropy density profiles 
through the EoS, and calculate the spacetime-rapidity 
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FIG. 6. (Color online) Centrality dependence of the initial 
entropy dSi/dris in ^ysNN = 2.76 TeV Pb+Pb collisions at 
the LHC with different values of Ksa.t, P and a. The curves 
are normalized such that such that the total entropy in 0 — 5 
% centrality class is one. 


density of the total entropy, dSi/drjs, as 


dS^ 

dijs 


= J dxdyTos{x,y,To)^, 


(39) 


where s is the local entropy density. In our case the initial 
velocity is zero and 7 = (l — = 1. Figurej^shows 

the normalized initial entropy as a function of centrality. 
The lines show calculations with different values of Ffsat, 
P and a. The actual entropy varies as the parameters are 
changed but to better compare the centrality dependence 
in different cases, we have normalized the results in this 
figure such that dSi/drjs = 1 in the 0 — 5 % centrality 
class in each case. As one can read from the figure, the 
centrality dependence changes only slightly with different 
values of iFgat and /?, but the width of the nucleon gluon 
distribution affects it much more. These extremes, i.e. 
a = 0.60 fm or 0.20 fm are, however, not supported by 
the HERA/ZEUS data. Note that, when coupled with 
viscous fluid dynamics, the ATgat = 0.45 and /3 = 0.8 
case corresponds to the case 77 /s = paramS in the data 
comparison in Sec. V, see Table |VI[ 

Figure shows the initial eccentricities £2 and £3 for 
the same cases as above. In addition we show the ec¬ 
centricity from the usual Glauber model initial state, 
i.e. a mixture of the eWN and eBC initial densities, 
e oc fpbin + (1 - f)Pwn, with / = 0.16. Similarly to 
the initial entropy case, there is practically no sensitiv¬ 
ity on ATsat and /3, but a strong sensitivity on the value 
of a. The pQCD + saturation initial conditions give 
values of £2 that are significantly larger than those of 
the Glauber model, but the £3 values are very similar in 
Glauber model and pQGD+saturation initial conditions 
with cr = 0.43 fm, i.e with the a value obtained from the 



FIG. 7. (Color online) Centrality dependence of the initial 
eccentricity in -^snn = 2.76 TeV Pb+Pb collisions at the 
LHC with different values of Aaat, P and cr (solid and dashed 
lines). The Glauber-model case is shown for comparison (dot¬ 
ted lines). 


HERA/ZEUS fit. 

Figure shows the entropy weighted average satura¬ 
tion scale Psat as a function of centrality in Pb+Pb colli¬ 
sions at the LHG, computed for the same values of ATgati 
P and cr as in the previous figures. Fig. & shows the 
same for Au+Au collisions at RHIG. Again, we see that 
the gluonic width a has the largest effect on the centrality 
dependence (compare the dashed lines) while P and ATgat 
affect more the normalization of Psat- The opposite sys- 
tematics in P and ATgat can be understood from Eq. (251 
at the naive scaling limit: Psat ~ {K /where the 
NLO/LO K-factor K of Eq. (27) increases with decreas¬ 
ing p. We also see that the average saturation scales 
remain above 1 GeV for a very wide range of centralities 
both at the LHG and RHIG. 

Figure]^ shows the fraction of the initial dSi/drjs from 
the regions of the transverse plane where Psat > 1 GeV, 
both in Pb+Pb collisions at the LHG, and Fig. Ep the 
same in Au+Au at RHIG, computed for the same values 
of Algat, P and a as above. Note again that the case 
with Algat = 0.45, P = 0.8 will correspond to the paramP 
case in the data comparison ahead in Sec. V. This figure, 
together with Fig. indicates that pQGD + saturation 
indeed gives the dominant part of the initial conditions 
over a sufficiently wide range of centralities both at the 
LHG and RHIG, and that the additional phenomenology 
at the low-density edges of the system does not play a 
major role. 

In Fig. we show the entropy production due to the 
viscous effects in fluid dynamics for Pb+Pb collisions at 
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FIG. 8. (Color online) Average psat as a function of centrality in ^snn = 2.76 TeV Pb+Pb collisions at the LHC (a), and in 
^snn = 200 GeV Au+Au collisions at RHIC (b) with different values of Asat, /3 and a. 



FIG. 9. (Color online) Fraction of dSi/drja from the region psat > 1 GeV as a function of centrality in ^snn = 2.76 TeV Pb+Pb 
collisions at the LHC (a), and in ^ysNN = 200 GeV Au+Au collisions at RHIC (b). 


the LHC, computed for the r]js = 0.2 and paramS cases 
(cf. Fig. 0. Since our starting time for the fluid dynamics 
is relatively small, tq = 0.2 fm, the entropy production 
becomes sensitive to the QGP viscosity: hence there is 
significantly more entropy produced for paramS in cen¬ 
tral collisions where the initial temperatures are highest. 
As we can see in the figure, the entropy production is 
rather significant but especially for the parametrizations 
where the QGP viscosity remains below that in paramS 
(and which will also reproduce the experimental data 
best) it can still be regarded as a correction. In practice, 
in order to get the same multiplicity, e.g. in the most 


central collisions, with all the different p/s parametriza¬ 
tions, Asat is adjusted for each p/s{T) separately. 


IV. FLOW COEFFICIENTS AND 
CORRELATIONS 

Before comparing our results with the LHC and RHIC 
measurements, let us recapitulate the definitions of the 
various flow coefflcients and correlations discussed in the 
next section. The azimuthal parts of the transverse mo¬ 
mentum spectra are, traditionally, decomposed into the 
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are replaced by 
dNch 


dr]psdp^d4> 
2 


A?7p 


E 


A?7| 


sinh 


ps 


dnh 


mT,i 


dN, 


dydp^dcj) ’ 


(44) 


where the sum is over all the charged hadrons, mT,i = 
a/ mf + p^, and rrii is the mass of the hadron i. 


A. Event-plane method 


In addition, it is also possible to define the so called 
event-plane flow coefficients as 


n„{EP}(pT) = (cos [n {(j) - 4'„{EP}))]),^, 


(45) 


where 


FIG. 10. (Color online) Entropy production as a function of 
the initial dSi/dpaS in -^/snn ~ 2.76 TeV Pb+Pb collisions 
at the LHC, computed for rj/s{T) from paramS in Fig.[^and 
for rj/s = 0.20. 


Fourier components and their phases or event-plane 
angles flin. For a single event these can be defined as 

(40) 

where the angular brackets (• • • )0 denote the average 

^ i^dydp^j) J dydp'^dcj) 

Similarly, the pT^-integrated flow coefflcients are deflned 
as 

n„(y)e“'""(^) = (e“^)^.p,, (42) 

where the average is defined as 

, , fdN\~^ . 9 dN , , 

<■■■>♦.^. = ( 1 ^) J («) 


^-nlEP} = 

1 (46) 

— atan 2 {{w cos{n(f>)), {w sm{n(j)))^^p^) , 

with w being a weight factor, e.g. w = px- The prob¬ 
lem with the event-plane method is that, although here 
it coincides with the previous definitions if ffnlEP} is 
defined appropriately ('I'„{EP} = if n; = 1), in the 
experiments there is a finite number of particles in single 
event, resulting in a finite resolution in determining the 
event-plane angle. The finite event-plane resolution in 
turn introduces the ambiguity to the relation between 
the underlying flow coefficients and the measured 
event-averaged event-plane coefficients (n„{EP})ev In 
the high-resolution limit (n„{EP})ev (vn)evj nnd in 
the low-resolution limit (n„{EP})ev —>■ (Vn)ev^- In the 
presence of the flow fluctuations, these two averages are 
in general different. Typically, the real events are some¬ 
where between these limits, and a consistent comparison 
to the data requires that the calculated events are ana¬ 
lyzed similarly to the experiments imi, and even then 
the exact experimental configuration, e.g. non-uniform 
acceptance, that deviates from a theoretical perfect de¬ 
tector, can introduce ambiguity to the results |I28| . 
Thus, in this work, we do not consider the event-plane 
flow coefficients but rely on those obtained from the cu- 
mulants discussed next. 


From here on we drop the y from the arguments, as we 
are using the boost-invariant approximation, where the 
flow coefficients do not depend on the rapidity. In prac¬ 
tice, the pT-integration is never over the full pT-range, 
but different experiments have different p^-cuts in their 
analyses — a fact to be taken into account in the calcu¬ 
lations as well. Also, in the case of unidentified charged 
hadrons the rapidity y cannot be measured, but the spec¬ 
tra are averaged over some pseudorapidity range A.r]ps 
symmetric around y = 0. In this case the spectra above 


B. Cumulants 

The ambiguity problem associated with the event- 
plane method can be resolved by using the n-particle 
cumulants. For example, the two-particle cumulant is 
defined as the correlation 

n„{2}2 = f 

N 2 J d(j)id(t>2 

( 47 ) 
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where dN 2 /d(j)id(j )2 is the two-particle spectrum (sup¬ 
pressing the possible rapidity and px dependence), which 
can in general be decomposed as a sum of a product of 
single-particle spectra and a “direct” two-particle corre¬ 
lation S 2 { 4 > 1 , 4 > 2 ), 


dN2 

d4>id4>2 


dN dN 
d(j)i d(j)2 




(48) 


The direct correlations can result e.g. from a p-meson de¬ 
caying into two pions, and these correlations are usually 
referred to as non-flow contributions. Using Eq. (42), the 
event-averaged two-particle cumulant can be written as 




{VI + S2) 


1/2 

ev 


\^n/ev 5 


(49) 


where the last equality follows in the absence of the non¬ 
flow contributions, i.e. assuming that all the azimuthal 
correlations are due to the collective flow only. It turns 
out that the 2-particle cumulant always results in 
regardless of the event-plane resolution |128| . therefore 
resolving the ambiguity in the event-plane method. 

In our calculations, we use the single-particle spectra 
directly, i.e. we are not considering individual particles. 
Therefore, in our calculations the event-plane resolution 
is in principle (up to the numerical accuracy) infinite, and 
we do not need the corrections due to the finite event- 
plane resolution. Furthermore, even though we compute 
the hadron decays, they are done at the level of single¬ 
particle spectra, and thus all the direct correlations (non¬ 
flow) are absent in our calculations. We also note that 
typically in the experimental analysis the non-flow cor¬ 
relations are suppressed by choosing e.g. pseudorapidity 
gaps between the pairs of particles in Eq. (47). 

For these reasons, for our purposes it is sufficient to 
define the cumulants directly through the flow-only limit. 
The pr-integrated 2-particle cumulant flow coefficients 
are then defined as 


Vn{‘^} — (Vn) 


= /„ 2 \ 1/2 
ev 5 


(50) 


where the Vn for a single event follows from Eq. (42) 


and the angular brackets denote the average over all the 
events in a given centrality class. Similarly, the event- 
averaged pT-integrated 4-particle cumulant flow coeffi¬ 
cients are defined as m 


U„{4} = - {vi)ev) 


1/4 


(51) 


In addition to the u„{2} and u„{4}, we also study the 
three-particle cumulant U4{3} measured by STAR |130| . 
defined as 


„,{ 3 ( ^ , 52 ) 

Originally, the higher-order cumulants were introduced 
to suppress the non-flow correlations IMl, but after the 
full realization of the importance of the event-by-event 


fluctuations m it has become clear that different cu¬ 
mulants do not only have different sensitivity to non¬ 
flow correlations, but also measure different moments of 
the underlying probability distributions of the flow coef¬ 
ficients. 


C. Event-plane correlations 

Different correlations between the flow coefficients and 
the event-plane angles give a rich variety of observables 
that can provide independent further constraints to the 
properties of the strongly interacting matter. In this pa¬ 
per, we consider also the correlations between the event- 
plane angles of the different harmonics. In principle, 
one could define the correlations between the angles di¬ 
rectly as (cos(fci'I'i -I- • • • -I- nkn'i>n))ev, with the an¬ 
gles defined according to Eq. (|4^, but as was noted in 
Ref. |128| . this leads to a similar ambiguity related to 
the event-plane resolution as for the event-plane u„{EP} 
discussed above. For this reason it was suggested that it 
is better to define the event-plane correlations as 


(cos(fci^'i H-hnA:„4'„))sp = 

COs(fci^i H-h nkn'^n))ev 

\/• • • {Vn'""^)ev 

where the k^s are integers with the property “ 

0. This definition is actually equal to the low resolu¬ 
tion limit of the (naive) definition above. These correla¬ 
tions were recently measured by the ATLAS Collabora¬ 
tion m, by using both definitions. 


V. RESULTS 

A. Multiplicities, pt— spectra and average pr 

Once we have fixed the coefficients of the non-linear 
terms in Eq. 0 from the kinetic theory calculations, 
and the width a = 0.43 fm of the gluonic T„ from the 
HERA data, we have essentially four free parameters 
{ATsat)/3) BJ/FS, ? 7 /s(T)} in our model. As shown in our 
previous studies |53l 155] , the parameters iLgat and /3 are 
strongly correlated and a continuum of equally well work¬ 
ing pairs can be found, however so that the experimen¬ 
tal data slightly favors larger values of jd. For simplic¬ 
ity, to reduce the number of free parameters, we fix here 
(3 = 0.8, and choose the BJ-case for the pre-thermal evo¬ 
lution discussed in Sec. |IIID| We then tune the remaining 
parameter so that the charged hadron multiplicity 
dNch/drjps matches the ALICE measurement in the most 
central Pb+Pb collisions, i.e., in the 0 — 5 % centrality 
class at the LHC. It should be emphasized that no fur¬ 
ther tuning is done for other centralities at the LHC, or 
for any of the RHIC results. 
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FIG. 11. (Color online) Centrality dependence of charged hadron multiplicities in ^sjvjv = 2.76 TeV Pb+Pb collisions at the 
LHC (panel (a)) and 200 GeV Au+Au collisions at RHIC (panel (b)), computed for the five p/s{T) parametrizations shown in 
Fig.Experimental data are from ALICE |133| . STAR |51| and PHENIX |134| . 


As discussed in Sec. |H C| we consider the five different 
p/s{T) parametrizations shown in Fig.lH The viscous en¬ 
tropy production, different for each p/^T) case, needs to 
be compensated by (iteratively) adjusting ATgat for each 
parametrization. The obtained values of ATgat aro shown 
for each p/s parametrization in Table VI The resulting 
centrality dependence of the charged particle multiplicity 
in Pb+Pb collisions at the LHC is shown in Fig. and 
compared with the ALICE measurements |133| . As can 
be seen from the figure, our calculation matches very well 
with the measured data, and in practice all the five p/s 
parametrizations give an equally good agreement. 


TABLE VI. The values of Ksat for different p/s parametriza¬ 
tions. 


p/s 

0.20 

paraml 

param2 

param3 

param4 

-^sat 

0.63 

0.50 

0.75 

0.45 

0.64 


Once the parameters are fixed at the LHC, the ^/s-, 
centrality- and also A-dependences follow from the calcu¬ 
lation. The comparison of the corresponding calculation 
for Au+Au collisions at the top energy of RHIC is com¬ 
pared to the PHENIX |I34| and STAR measurements 
in Fig. [TTJr. As can be seen from the figure, the agreement 
with the calculation and experimental data is again very 
good. We emphasize that here also the multiplicity in the 
most central collisions follows from the calculation, i.e., 
we do not change ATgat with ^/s or A. Because p/s{T) is 
also by definition independent of ^/s and A, we are now 
in principle equipped to predict the multiplicities in any 
other collision systems, provided that the fluid dynamics 
and pQCD + saturation pictures are valid. 

The comparison of the calculated pT-spectra of charged 
hadrons with the ALICE measurement |135| in Pb+Pb 
collisions at the LHC is shown in Fig. [T^ , and the 
corresponding comparison with the STAR |136] and 


PHENIX |137j data in Au+Au collisions at RHIC is 
shown in Fig. [T^d. As long as the multiplicities are well 
described, the pT-spectra are quite insensitive to the p/s 
parametrizations. In fact, the most important parame¬ 
ters that dictate the behavior of the pT-spectra are the 
kinetic and chemical freeze-out temperatures Tuec and 
Achem- While the multiplicity ratios of the identified 
hadrons, e.g. the pion-to-proton ratio, are best repro¬ 
duced with Achem ~ 150 MeV, it tends to give too fiat 
PT-spectra, especially in the law-px region, where fluid 
dynamics is expected to work best. This is the reason for 
our choice of a rather high Achem = 175 MeV. Although, 
the proton yields are somewhat overpredicted with this 
choice, we can, however, get a good description of the 
low-pr region of the charged hadron spectra, which we 
consider here more important than a detailed descrip¬ 
tion of the hadronic chemistry. Inclusion of bulk viscos¬ 
ity could help to improve the overall agreement with the 
data, see e.g. Ref. |5?| . 

Figure |13^ shows the average pj- for pions, kaons and 
protons compared to the ALICE measurements m- 
The pions are at low-pr the most abundant particles, 
and the very good agreement of our results with the data 
reflects the fact that the low-p^ region of the p^-spectra 
is well enough described. The same conclusion holds for 
the average pt in Au+Au collisions at RHIC, shown in 
Fig. [I^ against the PHENIX |139| data. While in both 
cases the average px of pions is well reproduced, espe¬ 
cially the centrality dependence of the proton (px) does 
not come out correctly. Whether this could be cured 
by a more detailed account of the chemical reactions 
in the hadron gas in the fluid-dynamical calculation, or 
whether a full microscopic treatment is needed, remains 
an open question. Overall, the agreement with the low- 
px charged hadron spectra, and the very good agreement 
with the pion average px gives us confidence that the px- 
integrated bulk observables for charged hadrons can be 
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FIG. 12. (Color online) Transverse momentum spectra of charged hadrons in y/SNN = 2.76 TeV Pb+Pb collisions at the LHC 
(panel (a)) and 200 GeV Au+Au collisions at RHIC (panel (b)), in the same centrality bins as in Fig. |11[ computed for the 
five rj/siT) parametrizations shown in Fig.Experimental data are from ALICE |135| . STAR |136| and PHENIX |137| . For 
visibility, the curves and the data points have been shifted by increasing powers of 10. 



FIG. 13. (Color online) Centrality dependence of the average pr for pions, kaons and protons in ^Jsnn = 2.76 TeV Pb+Pb colli¬ 
sions at the LHC (panel (a)) and 200 GeV Au+Au collisions at RHIC (panel (b)), computed for the five p/s(T) parametrizations 
shown in Fig. Experimental data are from ALICE |138| and PHENIX |139| . 
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well described within our framework. 


B. Flow coefficients 

The viscosity does affect the multiplicities through 
the viscous entropy production, but this effect gets here 
compensated by the re-tuning of iflsat for each 77 /s 
parametrization. Also, once the multiplicities are repro¬ 
duced, the details of the pT-spectra are quite insensitive 
to the values of Ty/s. Therefore, these quantities do not 
give a direct access to the determination of Ty/s from the 
experimental data. The most direct constraint to the vis¬ 
cosity of the strongly interacting matter comes from the 
azimuthal structure of the hadron spectra. 

The computed 2-particle cumulant u„{2} for charged 
hadrons at different centralities in Pb+Pb collisions at 
the LHC are shown in Fig. against the ALICE 
data |140| . For the definitions, see Sec. IV B As the fig¬ 
ure verifies, all the parametrizations of Fig. [^reproduce 
(by construction) the Vn’s at the LHC up to 40-50% cen¬ 
tralities very well. Thus, the LHC Vn data alone do not 
allow to distinguish between the different rj/s temper¬ 
ature dependencies to such a precision. More notable 
differences appear only in the more peripheral collisions, 
where the uncertainties related to the fluid dynamics and 
its applicability, as well as to the initial state calculation, 
are large. 

One can also note that the higher harmonics measured 
at the LHC do not give directly additional constraints to 
the temperature dependence of the viscosity. The ratio 
of V 3 or Vi to the elliptic flow coefficient V 2 , however, 
depends strongly on the initial conditions, through the 
ratio of the initial eccentricities e^lsn- Therefore, the 
higher harmonics give an indirect constrain to the ry/s, 
by restricting the possible initial states, see Ref. |148| . As 
seen in the figure, our approach with pQCD + saturation 
initial conditions describe the u„’s very well. 

So far, the u„’s at the LHC give at most the upper 
limit for the minimum oirj/s (corresponding to the con¬ 
stant 'qjs = 0 . 20 ), but even with these choices it varies 
between 7y/s|min = 0.08 and 0.20, with a possibility that 
even smaller Ty/smin could be tuned to fit the data. Fur¬ 
thermore, the location of the minimum is not constrained 
either. For the low- and high-temperature rj/s the uncer¬ 
tainties are even larger than for the minimum. It is then 
clear that further constraints are needed in order to pin 
down the temperature dependence of rj/s. 

A simultaneous analysis of other collision systems can 
provide further independent constraints for r]/s(T). As 
discussed in Refs. Ham], the viscous suppression of Vn ’s 
depends differently on the temperature dependence of 
r]/s{T) at different collision energies. In ^/s//// = 200 
GeV Au+Au collisions at RHIC the 7 ;„’s are practically 
independent of the high temperature, T 3> Tc, shear vis¬ 
cosity. At higher energies the high-temperature viscosity 
becomes gradually more important, while the influence 
of the hadronic viscosity decreases. 


In Fig. 14 3 we show the computed U2{2}, 773 ( 2 }, 


and 774 ( 3 } for charged hadrons in ^/snn = 200 GeV 
Au+Au collisions at RHIC compared to the STAR 
data |130L 114111142| . As one can read from the figure, 
the same r]/s{T) parametrizations that gave an equally 
good fit to the Vn data at the LHC are now clearly sep¬ 
arated, demonstrating that the simultaneous RHIC and 
LHC analysis of 77„’s can be used at least to rule out some 
temperature dependencies. Here, especially, the paramA 
with a large hadronic viscosity fails to describe the data. 
Overall, the best agreement with the data is obtained 
with a constant rj/s = 0.20 and rj/s from paraml with 
the minimum at T = 150 MeV. 


C. Flow fluctuations 

A proper event-by-event description of heavy-ion col¬ 
lisions collisions should not only reproduce the event- 
averaged 77„’s but also their EbyE probability distribu¬ 
tions P{vn). As we show here, and as earlier reported in 
Ref. [24], it turns out that the probability distributions 
of the scaled 77„, defined as 


5vn = 


77n (77n)e 

(77n)ev 


(54) 


are essentially independent of the details of the fluid dy¬ 
namical evolution, but depend only on the corresponding 
eccentricity fluctuations of the initial state. Therefore, 
the current LHC data on P{vn) provide a direct con¬ 
straint for the initial states m such as we compute 
here. 


In Figs. 15 1 and 15 3 we show the computed P{Sv 2 ) 
fluctuation spectra compared to the ATLAS data [55] in 
the 5—10 % and 35 — 40 % centrality classes, respectively. 
The pQCD+saturation initial state in this figure is com¬ 
puted with rj/s = 0.20. For comparison, we also show 
a calculation with the usual Glauber initial condition, 
where the energy density is proportional to a linear com¬ 
bination of a binary collision density pbin and a wounded 
nucleon density Pwn, i-e. e oc fphin + (1 - /)Pwn with 
/ = 0.15 and ry/s = 0.10 to approximately match the 
measured centrality dependence of the multiplicity and 
772 . The probability densities of the scaled eccentricities, 
P{Se 2 ) are also shown in the figures. 


As seen in the panel (a) of Fig. 15 in the near cen¬ 
tral collisions the scaled V 2 distribution follows closely 
the distribution of the scaled £2 with both the EKRT 
and Glauber initial states. The pQCD-based initial con¬ 
ditions give a very good description of the ATLAS data, 
while the Glauber initial conditions result in a too wide 
distribution. In mid-peripheral collisions, shown in the 
panel (b) of Fig. 15 the EKRT initial conditions give still 
a good description of the ATLAS data and the Glauber 
result is still too wide. However, as clearly seen in the 
figure, the scaled 77„ distributions do not anymore fol¬ 
low the eccentricity distribution, but the V 2 distributions 
are visibly wider than the £2 distributions, concretely 
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FIG. 14. (Color online) Centrality dependence of the flow coefficients ii„{2} from the charged hadron 2-particle cumulants in 
y'sjviv = 2.76 TeV Pb-t-Pb collisions at the LHC (panel (a)), and the coefficients V2{2}, V 3 { 2 }, and V 4 { 3 } from the charged hadron 
2- and 3-particle cumulants in 200 GeV Au-|-Au collisions at RHIC (panel (b)), computed for the five rj/s{T) parametrizations 
shown in Fig. Experimental data are from ALICE [140] and STAR [TMl ITiTl fT42| . 



FIG. 15. (Color online) Panel (a): Fluctuation spectra of the final-state V 2 of charged hadrons (solid curves) and of the initial 
state £2 (dashed) in the 5 — 10 % centrality class in ^snn = 2.76 TeV Pb-t-Pb collisions at the LHC, computed with the pQCD 
+ saturation initial states and rj/s = 0.20, and with the Glauber-model initial states using rj/s = 0.10. The experimental data 
are from ATLAS |28| . Panel (b): The same but for the 35-40% centrality class. 


demonstrating the necessity of fluid dynamics in describ¬ 
ing the detailed response to the initial eccentricities, see 
also Ref. |145] . The fluctuation spectra of the higher 
harmonics and V 4 are also well reproduced with the 
pQCD+saturation initial conditions, but they do not 
show similar sensitivity to the initial conditions as the 
V 2 fluctuations. 


Figure 16 shows the P{Sv 2 ) distribution of charged 
hadrons in the same 35—40 % centrality class with pQCD 
+ saturation initial conditions as panel (b) of Fig. 15 but 
with three different rj/s(T) parametrizations: rj!s = 0.20, 
77 /s = paramA, and p/s = 0. As can be seen from the fig¬ 
ure, the final Sv 2 distribution is the same with all three 
p/s parametrizations. This is true even in the perfect 
fluid limit p/s = 0. This shows that even if the fluid 


dynamical evolution plays a crucial role in getting the 
final V 2 distributions correctly reproduced in the periph¬ 
eral collisions, they are still a good probe of the initial 
conditions, because they do not depend on the details of 
the fluid dynamical evolution. 

Then, a very interesting question is how directly the 
final-state V 2 distribution can reflect the initial state £2 
distribution (and vice versa). If V 2 and £2 are, to a suf¬ 
ficient approximation, linearly correlated, V 2 oc £ 2 , then 
the scaled distributions P{Sv 2 ) and P{Se 2 ) are naturally 
identical. As seen from the panel (a) of Fig. 15 this is 
the case in central collisions. However, as noticed from 
the panel (b), the distributions are not anymore the same 
in peripheral collisions, indicating that there must be de¬ 
viations from the linear relation. What complicates the 
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FIG. 16. (Color online) Fluctuation spectrum of V 2 of charged 
hadrons in the 35 — 40 % centrality class in ^snn = 2.76 TeV 
Pb+Pb collisions at the LHC, computed with the pQCD + 
saturation initial states and with two different parametriza- 
tions of ri/s{T) and also using ideal fluid dynamics, rj/s = 0. 
The experimental data are from ATLAS |28| . 


initial state extraction from the V 2 fluctuation spectrum 
further is that the £„ = are not actually sufficient 
to determine the full angular structure of the initial den¬ 
sity profile, but in principle all of the em,n coefficients, 
defined in Eq. (351, are needed. 

In the left panels of Fig. we show the probability 
distributions of 6 v 2 , Se 2 , and fei ^2 in the 5 — 10 %, 35 — 40 
%, and 55 — 60 % centrality classes, obtained with the 
pQCD + saturation initial conditions and rj/s = 0.20. 
The middle panels show the correlation between V 2 and 
£ 2 ) and the right panels show the correlation between 
V 2 and £i^ 2 - In the 5 — 10 % centrality class all three 
distributions are practically the same. The linear relation 
between V 2 and £2 holds very well, thus the corresponding 
fluctuation spectra fall on top of each other. As the top 
right panel indicates, the correlation between V 2 and £ 1^2 
is visibly weaker, but the average V 2 computed at a fixed 
£ 1.2 still grows linearly with £ 12 , so that again P{Sv 2 ) ~ 
P((5£i,2). 

In the 35—40 % centrality class, the (u 2 , £ 2 )-correlation 
is still very strong, but there is already a clear deviation 
from a linear correlation, and as a result the V 2 and £2 
distributions are not anymore the same. However, the 
(u 2 ) £i, 2 )“Correlation is similar to the one in the near¬ 
central collisions, and the scaled V 2 distribution is prac¬ 
tically the same as the scaled £ 1^2 distribution. In even 
more peripheral collisions, i.e., in the 55 — 60 % cen¬ 
trality class, the (u 2 , £ 2 )~correlations show even stronger 
deviations from a linear correlation, and there is a slight 
deviation from the linear (u 2 , £i, 2 )^correlation as well. 

Overall, the (u2, £2)“Correlations are somewhat 
stronger than those of (u2,£i,2) but exhibit a strong 
non-linear behavior in more peripheral collisions. On 
the other hand, the (^2)£i,2)“Correlations stay more 
linear, and in central to mid-peripheral collisions the 
scaled V 2 distributions follow closely the scaled £1^2 


distributions, but in more peripheral collisions also 
they start to deviate from each other. Based on the 
middle and r.h.s. panels we can also deduce why 5 e 2 
spectrum in peripheral collisions becomes narrower than 
that of < 5 £i, 2 : for the averages (£ 2 ) > (£ 1 , 2 ) but the 
rare largest fluctuations are about the same magnitude, 
which for such largest absolute fluctuations means that 
£2 — (£ 2 ) < £ 1,2 — (£ 1 , 2 ), and for the scaled fluctuations 
even more strongly 5£2 < (£ 1 , 2 )- 

At the moment, it is not clear whether one could find 
a more specific definition of the eccentricity that would 
always be able to predict the V 2 distributions, or if the 
non-linear correlations remain inevitably a necessary part 
of the analysis. However, we emphasize that in a full 
fluid-dynamical calculation as presented here, the differ¬ 
ent definitions of the initial state eccentricities do not 
play a role in obtaining the final-state observables: The 
agreement between the ATLAS data and our calculations 
is very good, systematically over a wide range of central¬ 
ities. 


Another way to get an access to the flow fluctuations 
are the flow cumulants. As discussed in Sec. |IVBl if the 
non-flow contributions to the flow coefficients can be sup¬ 
pressed by the pseudorapidity gaps, the essential differ¬ 
ence between u„{2} and u„{4} is that they measure the 
different moments of the probability distribution P(u„). 
In principle, the full set of cumulants provides the same 
information as the probability distributions themselves. 
Thus, if we describe the U 2 { 2 } measurements simultane¬ 
ously with the full probability distributions, and the 
non-flow contributions are small, we should also agree 
with the U 2 { 4 } measurements. This turns out to be the 
case. 


In Fig. 18 1 we show the U 2 { 4 } of charged hadrons 
in Pb+Pb collisions at the LHC with different 77 /s 
parametrizations, against the ALICE data |I401II46] . For 


comparison we also show the V 2 { 2 } results from Fig. 14 


As one can see, the agreement with both measurements is 
equally good. Figure 18 3 shows the corresponding U 2 { 2 } 
and U 2 { 4 } in Au+Au collisions at RHIC compared to the 
STAR data |141] . The measurements of the full probabil¬ 
ity distributions are not currently available at RHIC en¬ 
ergies, but the fact that those rj/s parametrizations that 
give a good agreement with the U 2 { 2 } measurements also 
give an equally good agreement with U 2 { 4 } measurements 
already indicates that also at RHIC the main features of 
the probability distributions are correct in our approach 
with the pQCD + saturation initial conditions. 

If the flow fluctuations are approximately Gaussian, 
then U 2 { 4 } is approximately equivalent to the Uti{RP} 
determined with respect to the reaction plane (the calcu- 
lational plane whose x axis is along the impact parame¬ 
ter) |147| . In Fig. 18 1 we show also 7 ;„{RP} with 77 /s from 
paraml and paramA. As one can see, 7 ;„{RP} and u„{4} 
agree very well approximately up to the 40—50 % central¬ 
ities. Looking back at the left panels of Fig. [T^ this result 
is expected, since towards peripheral collisions the fluc¬ 
tuation spectrum exhibits more clearly a non-Gaussian 
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FIG. 17. (Color online) Left panels: Probability distributions of the charged hadron Sv 2 , and of the initial state Se 2 and & 12 , 
in the 5 — 10 % (top), 35 — 40 % (middle), and 55 — 60 % (bottom) centrality classes in ^/snn = 2.76 TeV Pb+Pb collisions 
at the LHC, computed with the pQCD + saturation initial states. The experimental data is from ATLAS |28| . Middle panels: 
The correlation between V 2 and £2 as a two-dimensional histogram. Right panels: The correlation between V 2 and £ 12 . The 
white lines in the middle and right panels are cubic polynomial fits, to guide the eye. The statistics for these figures was 15k 
events for each centrality class 


behavior. 


D. Event-plane correlations 

Because fluid dynamics is a non-linear theory, there is 
no reason to expect that the linear relation e.g. between 
the eccentricities and flow coefficients, oc e„, holds in 
general or even that is created by a non-linear response 
to the £„ alone. In reality, the different u„’s or l^n’s do 
not evolve independently, but are correlated with each 
other, e.g. a large V 2 can create a large U 4 even if the 
initial £4 is zero. The evidence for this can be clearly 
seen in the measured event-plane correlations |13211148] , 
which show a strong correlation between various event- 
plane angles iPn. 

Even though the correlation between the initial eccen¬ 
tricities creates correlations between u„’s through a lin¬ 
ear relation oc £„, even the signs of the measured 
correlations cannot be reproduced by this assumption. 


A generic behavior of the correlations can be explained 
by a linear response between the eccentricities defined 
through cumulants and u„’s, but quantitatively the 
magnitude of the correlations indicates that a non-linear 
fluid dynamical evolution is essential to reproduce the 
measurements, see Ref. m- Furthermore, and most 
importantly for the present study, the event-plane corre¬ 
lations give independent constraints to the initial state 
and transport coefficients, even if the viscosity is tuned 
to reproduce the V 2 data |151| . 


In Fig. we show various event-plane correlations in¬ 
volving two different event-plane angles T'n, defined by 
Eq. (531, in Pb+Pb collisions at the LHC, compared to 
the ATLAS measurements |132| . As can be seen from 
the figure, the different ij/s parametrizations that give 
an equivalent agreement with the data at the LHC, 
can be clearly distinguished by the correlations. Only 
two cases, ij/s = 0.20 and rj/s = paraml, give a good 
agreement with the ATLAS data. Only in the peripheral 
collisions (40 — 50 % centrality class and more peripheral) 
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FIG. 18. (Color online) 2- and 4-particle cumulant flow-coefficients, Vri{2} and of charged hadrons in ^snn = 2.76 

TeV Pb+Pb collisions at the LHC (a), and in 200 GeV Au+Au collisions at RHIC (b). The ?;„{4} results are divided by 2 for 
clarity. The dashed lines show the Vn calculated with respect to the reaction plane (RP). The data are from ALICE [140[ I146| 
and STAR |141| . and the corresponding pr ranges are indicated. 
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FIG. 19. (Color online) Correlations of two event-plane angles for charged particles in ^sjvat = 2.76 TeV Pb+Pb collisions at 
the LHC, compared with the ATLAS data |132| . 


the correlations involving T'g are not reproduced. A fur¬ 
ther discussion on how viscosity affects the correlations 
is given in the next section. 

The ATLAS Collaboration has also measured correla¬ 
tions involving three different event-plane angles m- 


As shown by Fig. these are equivalently well repro¬ 
duced in our framework by the same two parametriza- 
tions of T]/s as the two event-plane angle correlations 
above, but do not provide any further constraints to our 
setup so that rj/s = 0.20 and r]/s = paraml parametriza- 
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FIG. 20. (Color online) Correlations of three event-plane angles for charged particles in ^^snn = 2.76 TeV Pb+Pb collisions 
at the LHC, compared with the ATLAS data |132| . 


tions could be further separated by these measurements. 
It is to be emphasized that the same rj/s parametriza- 
tions that give the best fit to data at RHIC, also gives 
the best fit to the LHC event-angle correlators. 

In Au+Au collisions at RHIC the U4{3} measurement 
by the STAR Collaboration is actually similar to the 
event-plane angle measurement, as it involves also a cor¬ 
relation between the angles II '2 and II' 4 , see the defini¬ 
tion Eq. (52l. This particular measurement, shown in 
Fig. 14 D is also well described by the rj/s = 0.20 and 
77 /s = paraml parametrizations. 

Finally, we note that typically the required statistics 
(number of events) for the correlators is much higher 
than for the Vn coefficients themselves, and it also de¬ 
pends strongly on the strength of the correlation. For 
example, for the correlation between 472 and 473 , which 
is almost zero in Fig.[T9}I, the ATLAS Collaboration mea¬ 
sures clearly a positive value, while our current statistics 
(20fc events for each 77 /s parametrization) is not sufficient 
to accurately calculate such a small correlation but the 
statistical errors are larger than the signal itself. 


VI. DISCUSSION 

The dissipative suppression of the final azimuthal 
asymmetries of the spectra is a result of a combination of 
the dissipative effect into the fluid dynamical flow field. 


generated during the evolution, and the magnitude of the 
shear-stress tensor at the decoupling, i.e. the magnitude 
of the 6 f corrections to the equilibrium distributions in 
Eq. ([^. The relative contribution of these two effects 
depends on the 77 /s parametrization and collision energy. 

In order to illustrate the effects of Sf, we show in 
Fig.|2T}i the centrality dependence of V 2 at the LHC, cal¬ 
culated with three different constant 77 /s values ( 77 /s = 
0.1, 0.2 and 0.3), and with 77 /s from the parametrization 
paramA which has a large viscosity in the hadronic phase. 
The full results are shown with solid lines, and the results 
without the 5f contribution with dashed lines. Figure 
[ 21) 3 shows the same, but for 774 . Note that for these 
checks of the ( 77 , 5f) systematics we do not include the 
decay contributions, so the solid lines for 77 /s = 0.2 and 
paramA are not exactly the same as in Fig. but serve 
the purpose here. As one can see from the figures, the 
relative size of the Sf contribution increases with increas¬ 
ing rj/s, and also from central to peripheral collisions. In 
addition it is larger for higher harmonics, i.e. relatively 
larger for V 4 than for V 2 - 

In this work, we have tested different temperature- 
dependent parametrizations of 77 /s against the flow coef¬ 
ficient data from A+A collisions at RHIC and the LHC. 
First we noticed that the Vn measurements at the LHC 
alone do not give strong constraints on the temperature 
dependence oi rj/s but all our different parametrizations 
give an equally good agreement with the LHC data. We 
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FIG. 21. (Color online) Left: Centrality dependence of the flow coefficients V2{2} from the charged hadron 2-particle cumulants 
in ^/sn~n = 2.76 TeV Pb-|-Pb collisions at the LHC, for 4 different tj/s cases, with (solid curves) and without (dashed) the Sf 
corrections. Right: The same for Vi. The experimental data is from ALICE |14n| and the pr interval is indicated. 


emphasize that this is not trivially so, as the final az¬ 
imuthal asymmetry is generated in different ways with 
the different rj/s{T) parametrizations. This can be seen 
by comparing rj/s = 0.20 and rj/s = paramA curves in 
Fig. 1^ Both these parametrizations are tuned to re¬ 
produce the Vn data at the LHC, but the 6 f contribu¬ 
tion is significantly larger with rj/s = paramA due to the 
larger hadronic viscosity. Therefore, in order to repro¬ 
duce the data, the viscosity effects during the evolution 
need to be weaker in this parametrization compared to 
the ??/s = 0.20 case. It turns out, however, that even if 
the relative contribution from 5 f and from the evolution 
is different, the centrality dependence of is very closely 
the same with both parametrizations. One can see the 
differences only in very peripheral collisions where one 
has to be cautious about the applicability of fluid dy¬ 
namics. Therefore, the current Vn measurements at the 
LHC alone cannot reliably distinguish between the differ¬ 
ent temperature dependencies ofp/s at the level depicted 
in Fig. or in other words, they cannot be used to dis¬ 
tinguish the Sf contributions from the dissipative effects 
in the spacetime evolution of the flow field. 

For Vn both contributions, Sf and the dissipative ef¬ 
fects in the evolution, work in the same direction, i.e., 
both suppress the flow coefficients. Interestingly, the 
same is not true for the event-plane correlations. While 
Sf still suppresses the correlations, increasing the vis¬ 
cosity during the evolution can enhance the correlations. 
This is can be seen in Fig. where we show the event- 
plane correlations with the same 77/s parametrizations, 
again with and without Sf, as in the previous figure. In 
particular, one can see that the correlation between ’A12 
and 474 gets clearly stronger when 77/s is increased from 
0.1 to 0.3. The Sf contributions for these correlators 
remain small in the near-central and semi-peripheral col¬ 
lisions for all these 77/s parametrizations. Towards more 


peripheral collisions, however, the effect of Sf sets in, 
very quickly decorrelating the angles. 

Because of their different dependence on the viscosity, 
the event-plane correlations offer complementary infor¬ 
mation about the temperature dependence of 77/s. More¬ 
over, the weak dependence of (cos {N (472 — 474))) on the 
Sf in central and mid-peripheral collisions gives confi¬ 
dence that these correlations actually probe the dissi¬ 
pation during the evolution. Furthermore, the relative 
Sf contribution to also changes with collision energy, 
e.g., the Sf contribution is generally larger at RHIC 
energy. Therefore, it is remarkable that the same 77/s 
parametrizations that give the best agreement with the 
LHC correlation data also give the best agreement with 
the Vn data at RHIC. 

One should, however, keep in mind that large <5/ is a 
result of large values of inverse Reynolds number R~ ^ = 
TT^'^/Pg at the decoupling, which means that the system is 
not close to local thermal equilibrium, and the larger the 
the less reliable the fluid dynamical approximation 
becomes. Currently, it is not known to how large values 
of R~^ we can go in the current fluid-dynamical picture, 
so that we can still reliably calculate the evolution. 

Our results indicate that in order to keep the consis¬ 
tency with all the data shown here, the hadronic 77/s can¬ 
not be too large. At first this seems to be inconsistent 
with several microscopic calculations that show a strong 
increase of hadronic 77/s as temperature decreases, see 
e.g. Refs. HZllHSig. However, it should be noted that 
in our case, below the chemical freeze-out temperature 
Lchem = 175 MeV, the entropy density in 77/s is not an 
entropy density of the system in full chemical equilib¬ 
rium. Therefore, the comparison to microscopic calcula¬ 
tions, see e.g. Refs. |152H155] . which typically assume a 
full chemical equilibrium at all temperatures, would re¬ 
quire an estimate of how r]/s(T) = 77/spce(F) is related 
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FIG. 22. (Color online) Correlations of two event-plane angles for charged particles in = 2.76 TeV Pb+Pb collisions at 

the LHC, computed with (solid) and without (dashed) the Sf corrections. The experimental data is from ATLAS |132| . 


to the full equilibrium jy/scE- In order to estimate the 
magnitude of this difference, we show in Fig. 23 our rj/s 
parametrizations scaled with the ratio of entropy densi¬ 
ties in chemically frozen system and system in chemical 


equilibrium, i.e. rj/scE = (v/spce) x j • At least in 

a simplified hadron gas t] itself depends only weakly on 
the chemical composition |156| . and the main difference 
between jy/spcE and rj/scE is due to the change in the 
entropy density. The original parametrizations are shown 
as dashed curves. As one can see, the entropy densities 
of the two systems at low temperatures are significantly 
different. For example, the constant rj/s = 0.20 scaled by 
the entropy ratio, is very close to the original non-scaled 
77 /s = paramA parametrization, which is the one with 
the highest hadronic viscosity. 



T [MeV] 


VII. CONCLUSIONS 

In this paper, we have developed an event-by-event 
framework of the NLO-improved pQCD + saturation + 
viscous fluid dynamics model |55| . The main conclusions 
from the new EbyE EKRT framework are the following: 

1) We have now systematically tested the approach and 
successfully challenged it against a multitude of LHC and 
RHIC data. The centrality dependence of multiplicities, 
low-pT spectra, flow coefficients at the LHC and RHIC, 


FIG. 23. (Color online) Parametrizations of the tempera¬ 
ture dependence of the shear-viscosity to entropy ratio, scaled 
by the entropy density ratio of chemically frozen and chemi¬ 
cal equilibrium system. The dashed curves show the original 
parametrizations of Fig. 


and even the event-plane angle correlations at the LHC 
all come out in a beautiful agreement with experimental 
data. Especially the measured probability distributions 
of 6 v 2 at the LHC offer a stringent test for the computed 
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pQCD + saturation initial states. We have also demon¬ 
strated the necessity of fluid dynamical evolution in de¬ 
scribing the full centrality dependence of the measured V 2 
fluctuation spectra. The multiobservable analysis which 
is performed simultaneously for the LHC and RHIC, to¬ 
gether with our systematic fluid-dynamical cross-checks, 
suggests that the EbyE EKRT framework works remark¬ 
ably well for collisions up to 40... 50 % centralities. 

2) At the same time, as the main goal of this paper, we 
obtain improved constraints to the QCD matter r]/s{T). 
We tested several parametrizations of rj/s{T), all tuned 
to reproduce the u„{ 2 } in the mid-central collisions at 
the LHC. In practice, the centrality dependence of the 
Vn coefficients at the LHC alone do not give strong con¬ 
straints to the temperature dependence of r]/s, but all 
our parametrizations shown in Fig. give an equally 
good agreement with the data. The differences show 
only in peripheral collisions, where the uncertainties of 
the framework also grow large. A simultaneous analy¬ 
sis of the flow coefficients at RHIC gives more stringent 
constraints, and of the parametrizations considered here 
the constant 77 /s = 0.20 and 77/5 = paraml, with a small 
hadronic viscosity and minimum 77 /s at T = 150 MeV, 
give an overall best agreement with the flow coefficients 
at the LHC and RHIC. Especially 77 /s = paramA with 
the largest hadronic viscosity gives too strong a suppres¬ 
sion of the flow coefficients at RHIC. 

3) The event-plane angle correlations which have been 
measured at the LHC, and which are here shown to probe 
especially the viscous effects in the space-time evolution 
of the QCD matter, provide most useful additional and 
also rather stringent constraints for r]/s(T). Remarkably, 
again the same r]/s(T) which gives the best agreement 
with the RHIC Vn data, reproduces also the LHC event- 
plane angle correlations best. To put a real statistical 
error bar onto r]/s{T) requires a full global analysis of 
the LHC and RHIC heavy-ion bulk data, see e.g. Refs. 
|157H159] . This is clearly beyond the scope of our study 
here but we consider the present paper as an important 
step towards such an analysis. 

It is good to look back at the main uncertainties of the 
framework presented here. Our NLO calculation for the 
minijet Et is - as an IR/CL-safe calculation and with the 
given PDFs, pq, Ay and /3 - rigorous. The saturation as 
we consider it here, is a conjecture but clearly it captures 
quite correctly the dominant features in the initial mini¬ 
jet production, from which we then compute the initial 
energy densities and formation times locally in the trans¬ 
verse plane. Our handling of the pre-thermal evolution 
from the local formation times to the starting time of the 
fluid dynamical simulation could in principle be improved 
by giving the initial minijet energy densities to the fluid 
dynamics as source terms at the locally varying forma¬ 
tion times. Then, however, it is not clear whether the 


used fluid-dynamical picture is still valid as the density 
gradients and additional entropy generation at the ear¬ 
liest stages of evolution would become even larger than 
what they are in the present study with tq = 0.2 fm. Al¬ 
ternatively, one could develop an EbyE model also for the 
minijet production and feed the minijets obtained in each 
event into a parton cascade description such as DAMPS 
USD], and extract the initial conditions for fluid-dynamics 
(including also the possible initial transverse flow now as¬ 
sumed to be zero) at a later time. On the fluid-dynamics 
side, the largest uncertainties are related to the treatment 
of the late hadronic evolution, e.g. chemical and kinetic 
decouplings. This might improve if one couples the fluid 
dynamics with a hadron cascade in the hadronic phase 
at high enough temperature. Then, however, one type of 
model uncertainties are replaced with uncertainties re¬ 
lated to the matching conditions at the switching surface 
and uncertainties related to, e.g., the applicability of the 
cascade for very dense hadron systems, and also to the 
many unknown scattering cross sections one is forced to 
assume in such a simulation. 

The evident next step in our NLO-improved pQCD + 
saturation EbyE framework is to consider also the dy¬ 
namical fluctuations of initial gluon densities in the col¬ 
liding nuclei, which are then reflected as additional fluc¬ 
tuations of the saturation scale and hence of the com¬ 
puted initial energy densities. The inclusion of these 
fluctuations will improve our description of ultra-central 
heavy-ion collisions, and also allow us to study the ex¬ 
tremely interesting question of collectivity and flow p+Pb 
collisions at the LHC, see e.g. Refs. |16im65] . An inter¬ 
esting further question is the rapidity dependence of all 
the observables studied here. For this, one needs to de¬ 
velop a more complete EbyE framework by introducing 
a pQCD minijet event generator which is coupled with 
the determination of saturation in each event and which 
by construction also accounts for the different types of 
fluctuations. 
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